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INTRODUCTION 


The Stirling engine is currently the subject of intensive study by the 
Department of Energy. This engine has the potential to provide clean 
and efficient propulsive power for the automobile and for other high- 
way vehicles. In order to do so, however, several technical problems 
associated with such use of the Stirling engine must be solved. 

One of the technical problems that must be addressed is the rod seal. 
This seal has the function of separating the high pressure gas from 
the low pressure crank case oil. The seal can be either a sliding 
seal or a roll sock seal. The sliding seal is simpler, less expensive, 
and possesses a more gradual failure mode than the roll sock seal. 
However, the sliding seal, unlike the roll sock seal, is not hermetic. 
Additionally, the sliding seal is subject to wear and can produce un- 
desirable levels of friction force. 

The performance of the sliding seal with respect to leakage, wear, 
and friction can be improved from its present level. To do so re- 
quires that the details of the behavior of the seal in the Stirling 
engine be understood. Additionally, it is necessary to be able to 
treat the behavior of the seal in a quantitative manner. It is the 
objective of the work under the current NASA Contract DEN3— 22 to pro- 
gress substantially toward these ends. Specifically, the contract is 
directed at applying hydrodynamic and elastohydrodynamic theory to 
the rod seal. The contract is also concerned with the experimental 
determination of film thickness, fluid leakage, and power loss. The 
contract entails both producing correlation of experimental and 
theoretical results and developing tools appropriate for evaluation 
of rod seal behavior. 

The present document is an interim report of work completed in the 
first year of the two year effort. During the year tools for evalua- 
ting rod seal behavior have been developed. Analytically, a computer 
model of the elastohydrodynamic behavior of the rod seal for recipro- 
cating motion is available. The model permits the large initial 
radial squeeze, which is typically used for elastomeric seals, to 
exist. The model determines both the pressure distribution and the 
oil film thickness distribution in the seal/cylinder contact zone. 

In this, the calculations are made with time as the independent param- 
eter; consequently, the computer model provides film thickness and 
pressure results as a function of cyclic position. 


Experimentally, a rugged apparatus capable of providing film thickness, 
friction force, and seal leakage data has been completed. The apparatus 
contains a moving transparent cylinder (guided by precision hydrostatic 
bearings) and a stationary elastomeric test seal. A pressure gradient 
of 690 kPa (100 psi) can be applied across the seal. Frequencies from 
10 Hz to 50 Hz with a 0.254 m (1.000 inch) total stroke can be 
employed. Film thickness is measured with interferometry, 
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fluid leakage by level and pressure changes, and power (friction) loss 
by force cells. 

This report describes both the experimental and the analytical work 
completed to date. The experimental effort is discussed in the next 
section. Included in the section are descriptions of the experimental 
apparatus, of its capabilities, and of preliminary results from its 
use. Also included in the section is a discussion of the technique, 
developed under the contract work, for producing interferometry fringes 
at the seal/cylinder interface. The special attention given in this 
report to interferometry is due to the importance of film thickness 
measurement in the contract work. 

The analytical effort is covered both in the section following that 
for the experimental work and in the Appendices. The section presents 
the rod seal analysis which is implemented in the computer program 
given in the Appendixes. (A description of how the program is used is 
also given in the Appendixes.) In addition, the section presents results 
obtained from use of the computer program. 

The Conclusions and Recommendations section concerns the application 
during the ensuing year of the analytical and experimental tools 
developed during the past year and described in this report. 
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EXPERIMENTAL REPORT 


This section discusses the experimental apparatus that was designed 
and constructed during the past year. The capabilities of the device 
and the results obtained from its preliminary use are described. In 
addition, measurement of film thickness, accomplished by interferometric 
techniques, is considered. The separate treatment is given to inter- 
ferometry because of the importance of film thickness measurement to 
the contract work. 

The experimental device has not as yet been used to evaluate the be- 
havior of the rod seal or to produce experimental/ analytical correla- 
tions. Work in these areas is planned for the ensuing year and is 
described in the Conclusions and Recommendations section. 

The Experimental Apparatus 

The experimental apparatus is designed to facilitate the measurement 
of film thickness, seal friction force, and fluid leakage. The 
greatest experimental difficulties are associated with determining 
film thickness. Consequently, the design of the apparatus is based 
on the technique selected for this measurement. 

A. review of the available techniques for film thickness measurements 
was made in order to select that method most appropriate to the elasto- 
meric seal/cylinder wall interface. Some details of that review are 
presented in the subsection below entitled "Interferometry . 1 The 
result of the review was the selection of optical interferometry. 
Accordingly, the experimental apparatus contains a transparent cylinder. 
In addition, in order to minimize the difficulties of using the inter- 
ferometric technique, it is this cylinder, rather than the rod seal, 
that reciprocates. The use of a reciprocating transparent cylinder 
and a stationary elastomeric rod seal was the most important feature 
that affected the design of the device. 

Other features also affected the design. These features can be grouped 
in three categories: Measurement capabilities, performance capabilities 

and characteristics, and ease of use. Each of these categories is con- 
sidered below. 

In the measurement category, it is desirable to be able to determine 
the instantaneous seal friction force and the instantaneous pressure 
on the gas side of the seal. It is also desirable that some indica- 
tion be obtained of the temperature near the seal/cylinder interface. 
Finally, it is desirable for the apparatus to produce direct macro- 
scopic indications of oil and gas leakage. 

In the category of performance capabilities and characteristics, a 
rather lenghtly list of desirable features can be developed. Some 
of the more significant ones are: 
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. capability for a reasonable magnitude of pressure drop across 
the seal 

. compression (rather than tension) in the transparent cylinder 
(to utilize its strength characteristics) 

. minimum variation of the gas pressure during each cycle of the 
cylinder 

. precise guidance of the cylinder and minimum cylinder runout 
. infinitely variable speed in the range from 10 to 50 Hz 
. provision for water cooling 

. accurate alignment of the seal with respect to the transparent 
cylinder 

. long life span (low stress levels) 

. rigid mounting structure 
. few static seals 

. oil above the test seal (to insure fully flooded conditions) 

. vertical rod axis (to insure axisymmetric seal behavior) 

. balanced dynamic forces, and 

. synchronization of optics with cylinder motion. 

In the category of ease of use, significant features are: 

. ease of testing various types of seals 
. ease of replacing the seal and maintaining alignment 
. ease of removing the test seal from the transparent cylinder, 
and 

. ease of parts replacement. 

The experimental apparatus constructed during the past year was de- 
signed to incorporate the features in the three categories above. 

The design is centered around the transparent cylinder which recipro- 
cates with a 25.4 mm (1.000 inch) peak to peak displacement. The 
transparent cylinder has a 76 mm (3 inch) I. D. This I. D. 
was chosen large enough to prevent excessive difficulties in producing 
optical interferometry measurements, and small enough to keep the 
size of the apparatus (and the associated internal forces) reasonable. 
In addition, the cylinder bore was chosen such that the seal size is 
standard and relatively common. 

The overall apparatus is shown in the photograph of Figure 1. In the 
photograph can be seen the stationary plunger which holds the test 
seal and the oscillating assembly which contains the transparent 
cylinder. The photograph also shows, at the lower left, a pump and 
pressure tank. This pump and pressure tank provide high pressure oil 
to precision hydrostatic bearings which guide the reciprocating parts. 
The large motor in the center of the photograph is the prime mover 
which drives the apparatus. Its output speed can be varied contin- 
uously from less than 10 Hz to more than 50 Hz. The lattice cage 
which encloses the rotating and reciprocating parts provides a measure 
of operator safety. The framework for the apparatus is constructed 
of thick steel plate and I beams with relatively large cross sections. 
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The frame stands about 1.8 m (6 feet) tall. As pictured in 
Figure 1, the device has a mass of approximately 500 kg 
pounds mass) . 

Removal of the guard cage allows components of the drive system to be 
seen. This has been done in Figure 2, which shows the flywheel, two 
connecting rods, and 4 hydrostatic bearings with their associated 
plumbing. These components can also be identified by referring to 
Figure 3, which is a schematic drawing of the apparatus. 

The lower connecting rod is attached to a shaft which moves downward 
when the upper assembly (which contains the transparent cylinder) 
moves upward. The sole function of this lower shaft is to balance . 
the vertical forces in the rig. Each shaft is guided by two precision 
hydrostatic oil bearings which are of a four pad design and contain 
rubber seals at top and bottom. The shafts are constructed of alumi- 
num. They are attached to their respective connecting rods by needle 

bearings. 

Figure 4 shows some details of the upper end of the apparatus . A. 
schematic drawing which corresponds generally to this photograph is 
presented in Figure 5. The schematic shows both the exterior outline 
of this end of the apparatus (on the left in the figure) and some in- 
ternal details (on the right in the figure). 

Referring to both Figures 4 and 5, the hydrostatic bearings and the 
plunger which holds the test seal are mounted to a backing plate. 

This backing plate and the one for the lower shaft were assembled to 
the framework (against the machined tabs, Figure 4) and the bores of 
all four hydrostatic bearings were machined together. This process 
insured that the axes of all four hydrostatic bearings became co- 
incident upon assembly of the apparatus. During this machining pro- 
cess, the stationary plate was milled. This procedure resulted in 
the top surface of the stationary plate being perpendicular to the 
axis of the four hydrostatic bearings. Doweling of the stationary 
plate to the backing plate and of the hydrostatic bearings to the 
backing plate allowed for disassembly of the apparatus subsequent to 
the machining operation for insertion of the moving shaft. 

Several additional components in Figure 4 are worthy of discussion. 

The vertical pipe at the lower left of the photograph is the manifol 
for the hydrostatic oil supply. The three small pipes just to the 
left of center at the top of the photograph are the supply lines to 
the oil side of the seal, to the gas side of the seal, and for the 
water cooling. The return water leaves the apparatus through the 
threaded pipe at the top center of the photograph. 

The plunger is connected to the stationary plate with four bolts 
(three of which can be seen in Figure 4) . This arrangement can be 
readily described by referring to Figures 6 and 7. Figure 6 shows the 
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Figure 2 Overall View of Experimental Apparatus 
With Guard Cage Removed 
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apparatus with the plunger removed. The stationary plate is * 
top of the picture, the moving transparent cylinder (and its housing) 
are at the center of the picture, and the uppermost hydrostatic . 
bearing is at the bottom of the picture. Figure 7 is a view 

with the plunger installed. In Figure 7, the bottom of the stationary 
plate can be seen at the very top of this photograph. 

Figures 6 and 7 (as well as the schematic of Figure 5) indicate how 
the transparent cylinder is supported in its housing. The cylinder 
is inside the webs of this aluminum housing and is compressed axially 
by the 8 springs. Alignment is produced by boring the cylinder at 
the same time as the bottom face is cut. Not shown m the photograp s 
(but included in Figure 5) is an outer transparent cylinder. This . 
outer cylinder is also loaded axially by the 8 springs and can sustain 
690 kPa (100 psi) of gas pressure. Use of the outer cylinder resu ts 
in compressive, rather than tensile, gas loading on the inner 
cylinder. The apparatus can be run without the outer cylinder (for 
zero pressure drop across the seal) or with the outer cylinder ( or 
pressure drops up to 690 kPa (100 psi)). across the seal. Also not 
shown in the photographs (but included in Figure 5) is the containment 
for the gas volume. This containment is effected by attaching a 
Bellofram rollsock seal to the inside of the moving upper shaft. 

The other end of this rollsock seal is attached to a stationary plug 
which is in turn, attached to the backing plate. The design o 
the rollsock seal is such that the gas volume changes little during 
the reciprocating shaft motion. When the gas volume is connecte 
to an external gas reservoir, the peak pressure variations are esti- 
mated to be less than 1.70 kPa (0.25 psi) at 690 kPa (100 psi) gas 
pressure. 


Figure 7 shows several significant parts of the stationary plunger. 
Through the transparent Lexan cylinder, one can see, proceeding rom 
top to bottom, the felt wiper which contains the seal oil, plugged 
holes which are associated with the water cooling system, a volume of 
test oil (upper dark band), and the test seal (lower dark band). 

The parts can be seen more easily in Figure 8, which shows only the 

plunger . 

In Figure 8, the test seal can be seen near the lower right end of 
the plunger. This end of the plunger (the test seal holder) can be 
disassembled from the body of the plunger by removing the 4 bolts at 
the lower right of the picture. This is accomplished by firs 
straightening the locking tabs. Having the test seal holder replace- 
able allows various types of seals to be tested. For each ype o 
seal, a holder is required in which a groove appropriate to the seal 
has been cut. For the holder in Figure 8, the groove has been cut for 
the T seal shown in the photograph. 
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The plunger mounts to the stationary plate at the lower surface of the 
plate having the holes (Figure 8) . Just below this plate is a force 
cell (Kiag type 9071) for measuring the friction force produced by 
the seal. The connector for the lead is the white protrusion just 
above the cooling water supply hole. The force cell is preloaded 
between the plate and the remainder of the plunger body by a hollow 
tie bolt (the nut that controls the magnitude of the preload is at the 
upper left in Figure 8) . The assembly of the plunger to produce the 
preload is a one-time operation — the tie bolt portion of the plunger 
need never be disassembled at any time during all planned use of the 
apparatus. 

Use of the tie bolt to preload the force cell is desirable for several 
reasons. The first reason is that the hollow tie bolt is a convenient 
way to remove the cooling water from the plunger. The second reason 
is that with a properly designed tie bolt only one force cell need 
be used to measure seal friction. To see this, the following diagram, 
which shows the axial elastic behavior of the plunger, can be used: 



The diagram indicates four stiffnesses, those of the tie bolt, of the 
plunger below the plate, of the plunger above the plate, and of the 
force cell. A friction force, F, at the seal/cylinder interface must 
produce an axial displacement, d, of the test seal holder such that 


F 


^ k pa + k f k pb i # d 
k pa + \ V + k f j 
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Also, the force, F^, measured by the force cell is 


so that 


or 


F = 
c 


k k 

f P b .d 

k P b + k f 


k f k pb /(k pb + k £ > 


•s k pa /(k pa + V + k f k pb /<k pb + k f ) 


F " “b k pa (k pb + k f )/(k f k pb (k pa + V> + 1 

Now, if the design of the plunger is such that is small compared 

to k , the above expression reduces to 
pa’ 


F 

c 

F 


k b /k f + V k pb + 1 


which further reduces to F c /F = 1 if 1^/^ and l^/k are small com- 
pared to 1. The plunger was designed to meet these conditions so that 
the applied and the measured forces are nearly identical. Evidence 
that this design objective has been met is presented in the prelimi- 
nary results section below. 

Figure 8 can be used to point out other design features of the plunger. 
One such feature is the reduced diameter region below the test seal. 

This has the function of easing the installation of the test seal on 
the test seal holder. Another feature is the pressure transducer for 
measuring the small cyclic variations of gas pressure. This measure- 
ment can be made on a real time basis and can be used to subtract the 
gas pressure force variations from the force measured by the load 
cell. The lead from the pressure transducer can be seen (white wire) 
emerging from the plunger at the upper center of the photograph. The 
other two wires are thermocouple leads for measuring temperatures in 
the oil above the seal and in the test seal holder near the seal groove. 


The present test seal holder (that shown in Figures 1, 2, 4, 7 and 8) 
was designed to accept a T seal. This T seal consists of an elastomeric 
portion, one part of which contacts the cylinder, and two fiber back- 
up rings. The mounted seal is shown schematically in Figure 9. 
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Figure 9 Schematic of "T" Seal as Installed 
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Figure 10 is a photograph of the seal prior to installation on the 
test seal holder. The elastomeric material of this seal is Viton, a 
f luroelastomer . 

Preliminary Results 

In this section, results obtained from preliminary use of the experi- 
mental apparatus are presented. These results concern the operation 
of the apparatus rather than the behavior of the seal. The results 
show that the device performs well mechanically and that its force 
measurement capability is as intended. 

Upon completing the assembly of the rig, the first measurement made 
was of the stiffness of the hydrostatic bearings. A high degree of 
stiffness is desirable in order that dynamic forces not result in 
appreciable horizontal motion of the transparent cylinder with respect 
to the seal. The measurement was made by simply applying a steady 
horizontal force between the two hydrostatic bearings of the upper 
shaft. This force was measured with a spring scale. The displace- 
ment was measured with a dial indicator. The results showed the both 
upper hydrostatic bearings together have a horizontal stiffness of 

1.2 x 10^ N/mm (700,000 pounds per inch). The stiffness of 

each bearing is, therefore, at least 0.6 x 10^ N/mm (350,000 
pounds per inch) . 

To show what the effects of this stiffness are on the seal/cylinder 
interface, two rough calculations can be made. For these calculations, 
the following conservative assumptions are used: 

- The horizontal stiffness of each upper hydrostatic bearing is 

0.6 x 10^ N/mm (350,000 pounds per inch). 

- The horizontal dynamic force in the apparatus is 290 N 

(65 pounds). This force is calculated using the geometry of 
the connecting rod (the ratio of crank length to connecting 
rod length is 0.05), the maximum frequency (50 Hz), and the 
mass of the upper reciprocating shaft. This mass is less than 
4.5 kg (10 pounds mass). 

- The horizontal dynamic force is applied 0.13 m (5.0 inches) 
below the lower hydrostatic bearing of the upper shaft and these 
hydrostatic bearings are 0.21 m (8.1 inches) apart. 

- The seal is 0.22 m (8.5 inches) above the upper hydrostatic 
bearing of the upper shaft. 

The displacement of the transparent cylinder results from two effects — 
translation of the upper shaft due to the 290 N (65 pound) force, 
and rocking in the vertical plane due to the same 290 N (65 pound) 
force. The translation, t, can be computed simply from 
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t — =* 2.4 x 10 ^ mm (0.1 mils) 

2(0.6 x 10 8 ) 

This motion Is clearly negligible when compared to a typical seal 
preload. The displacement, r, due to rocking can be computed from 


r = 


290 


0.6 x 10 


8 


To. 13 + 0.105*1 |o. 22 + 0.105 1 
[ oTnT ][ 0.105 J 


_ 2 

= 0.17 X 10 Iran (0.67 mils) 

This motion is also negligible when compared to a typical seal preload. 

The next measurements to be made after assembly of the experimental 
apparatus were concerned with the smoothness of the device during 
operation. To obtain an evaluation of this, an accelerometer was 
placed on the backing plate at a height approximately corresponding 
to that of the seal. The sensitive axis of the accelerometer was 
oriented first in the vertical direction and then in the horizontal 
direction (such that the axis was in the plane of the backing plate) . 

The drive motor was operated at about 900 rpm. The plunger was removed. 

The results obtained are shown in Figures 11 and 12, where Figure 11 is 
for the horizontal acceleration and Figure 12 is for the vertical 
acceleration. Each figure gives an average spectrum made from 8 
spectrum analyses of the measured acceleration. The analyses were 
made in the 0-200 Hz range and therefore cover the flywheel frequency 
of 15 Hz as well as the drive motor frequency of about 30 Hz. 

The plots show that the vibration levels in the apparatus are relatively 
low. Reference to typical charts* used in vibration analysis indicates 
that both the horizontal and vertical vibration levels fall in the 
"smooth" range for machine operation. This shows that a satisfactory 
dynamic environment exists in the experimental apparatus for conducting 
the optical interferometry studies. 

The plots show that the horizontal acceleration is greater than the 
vertical acceleration, particularly in the regions near 60 Hz and 
30 Hz. This may have been anticipated, since the vertical, but not 
the horizontal, forces are balanced by the lower reciprocating shaft. 
However, these two frequencies are not directly associated with the 
rotation speed of the flywheel, which is about 15 Hz. 


*0ne such chart is the "Vibration Computer," a paper slide rule, pub- 
lished by Endevco Corp. 
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The origin of the acceleration signals at these frequencies is of 
interest, so additional spectra, Figures 13 and 14, were produced. 

For these spectra, all conditions were identical to those of Figures 11 
and 12; however, the coupling between the drive motor and the flywheel 
was disconnected. Consequently, the spectra of Figures 13 and 14 
were obtained under the condition of a stationary flywheel (and 
stationary connecting rods) but a rotating drive motor whose output 
speed was about 900 rpm (15 Hz) . 

Comparison of Figures 11 and 13 shows that most of the horizontal 
acceleration at 30 Hz, 60 Hz, and 180 Hz is not due to the flywheel 
and reciprocating masses; in fact, the only significant flywheel/ 
reciprocating mass signal is the one at about 53 Hz. Comparison of 
Figures 12 and 14 gives a similar result. For vertical vibrations, 
it can be observed that the flywheel/ reciprocating mass portion of 
the apparatus contributes essentially nothing to the vibration 
spectrum. 

One further plot was made. The intent of this plot was to determine 
whether the 30, 60, and 180 Hz signals were due to the drive motor or, 
perhaps, due to some electrical source. The plot, Figure 15, is an 
instantaneous spectrum made from the accelerometer with the drive 
motor switched off. The plot shows that significant portions of the 
60 Hz and 180 Hz signals are still present. These residual portions 
are, therefore, not associated with mechanical vibrations produced 
either by the drive motor or by the flywheel/reciprocating mass com- 
bination. This indicates that the drive motor contributes much less 
vibration to the experimental apparatus than Figures 13 and 14 suggest. 
Consequently, of the vibration levels of Figure 11, half of the 60 Hz 
and almost all of the 180 Hz signals are not mechanically produced. 

Also, most of the remainder of the 60 Hz signal and most of the 30 Hz 
signal are produced by the drive motor. 

One other mechanical aspect of rig operation is of concern. That 
aspect is whether a vertical force (i.e., the seal friction force) 
applied to the plunger at the seal is accurately measured by the force 
transducer. (The stiffness conditions for this to occur are consid- 
ered in the section above.) To evaluate the performance of the 
force transducer, a reference force transducer was attached to the 
bottom of the plunger at the hole on the vertical centerline. To 
this reference transducer was fastened a spring mass combination. 
Schematically, the arrangement is as shown in Figure 16. The spring 
mass combination was arranged so that the natural frequency was 
about 12 Hz. 

Displacement and subsequent release of the mass resulted in an 
oscillation which applied a sinusoidal vertical force to the reference 
cell. The force was also transmitted to the bottom of the plunger. 

Both the force measured by the reference force cell and that measured 
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Figure 16 Location of Reference Force Ceil 
and Spring-Mass with Respect 
to the Plunger 


by the force cell in the plunger were recorded by photographing the 
voltage outputs from their electronic drivers. These drivers were 
set so that each load cell had the same sensitivity 2.25 mV/N 
(10 mV/pound). The results are given in the photographs of Figure 17. 

The three sets of time traces in Figure 17 are for segments from three 
separate releases of the mass. The first set corresponds to a vibra- 
tion amplitude such that no spring bottoming occurred. The second set 
involved some bottoming near the start of the trace shown. The third 
set was made under a condition of relatively hard spring bottoming at 
the start of the oscillation. In each set of traces, two curves are 
given — the top one is for the reference force transducer and the 
bottom one is for the force cell in the plunger. The vertical and 
horizontal scales are as shown in the figure. 


All three photographs indicate very good correspondence between the 
forces measured by the two cells. The time traces are in phase, in- 
dicating that mass effects are not significant, even for the higher, 
frequencies that occurred. The magnitudes correspond well, suggesting 
that the required stiffness relationships for the tie bolt, plunger 
assembly and plunger force cell have been met. It therefore, appears 
that the force cell in the plunger does satisfactorily measure vertical 
forces (i.e., seal friction forces) that are applied to the test seal 
holder. 

Interferometry 

The accurate measurement of the film thickness at the seal/cylinder 
interface is central to the success of the current program. Film 
thickness measurements are the most sensitive means for performing 
analytical and experimental correlations. In addition, knowledge of 
film thickness can be used in studying seal friction behavior and 
can provide insight into how leakage can be controlled. 


Historically, many techniques of measurement have been used to ob- 
tain EHD film thickness. Successful methods of film thickness de- 
tection have been developed with the aid of optical equipment, capac- 
itive probes, inductive probes (eddy current devices), and resistance 
measurement sensors. The optical approach, however, is the only 
direct scheme which can provide a full point- by- point topographical 
film measurement within a reciprocating seal interface. In addition, 
the optical approach allows gradients within the contact region to 
be observed and film thicknesses as small as 2 microinches to be 
measured reliably. For these reasons, an optical approach to film 
thickness measurement has been taken and the moving transparent 
cylinder concept has been used. 

Application of optical procedures to measuring thicknesses such as 
those in elastohydrodynamic contacts invariably involves the inter- 
ferometry technique. This technique, which relies on the interference 
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between two coherent beams of light, is described beiefly below. 

Figure 18 shows a greatly expanded drawing of the seal/cylinder 
contact zone. The orientation of the seal and cylinder in the draw- 
ing is the same as that in the experimental apparatus — the section is 
made by passing a vertical plane through the vertical axis of the 
transparent cylinder. The drawing includes a light source which 
sends a light beam. A, toward the outer surface of the cylinder. This 
beam, which also lies in this vertical plane, has an intensity 1^. 

At the outer surface of the cylinder, this beam splits into two parts; 
one part, B, having intensity 1^, is reflected. The other part 

enters the cylinder wall and proceeds to the inner surface of the 
cylinder. At this surface, the beam splits such that one part is re- 
flected internally in the glass and the other part passes into the 
oil, to the surface of the elastomer, and then back to the inner wall 
of the glass. The viewer, looking in the direction indicated sees 
beams B, C, and D, where C is form the glass/oil interface and D is 
from the seal/oil interface. If certain physical conditions are met, 
interference between beams C and D will occur . One of these con- 
ditions is that the ratio of the thickness h to the wavelength of the 
incident light is h, 1, 1*2, etc. If the film thickness changes with 
axial position, the viewer will see alternating light and dark bands. 
These bands are termed interference fringes. It is evident that the 
axial spacing of these fringes gives directly a map of the changes in 
film thickness throughout the contact region. 

Several physical conditions must be met, in addition to the proper 
relationship between wavelength and film thickness, for observable 
fringes to exist. One of these conditions is that the beam inten- 
sities I and I. be approximately equal. Another condition is that 
c d 

the intensities I and I, be sufficiently high to allow them to be 
c d 

seen and photographed. A third condition is that 1^ be sufficiently 

small (or be somehow eliminated) so that it does not overwhelm the 
interference between C and D. Finally, the interface between the 
glass and the oil and the interface between the oil and the seal must 
be sufficiently specular that beams C and D are well defined. 

Satisfying the above conditions for an elastomeric seal is not straight- 
forward. The most difficult problem to overcome is the creation of 
a specular (optically smooth) surface at the oil/seal interface. 
Additional, but more readily overcome, problems are to minimize the 
effects of beam B, and to produce sufficiently intense beams C and D. 
Because of these problems, a little optical interferometry is done 
using elastomers such as rubber, and few reports of such work can be 
found in the literature. In those reports that exist (see, for example, 
[1] - [3]*), the problems are dealt with as follows: 


*Numbers in brackets denote references. 
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• Creation of a specular surface is accomplished by molding 
the elastomer in an optically smooth mold, thereby creating 
an optically smooth surface on the elastomer. 

. Creation of a specular surface by coating the elastomer with 
a thin specular material such as a metal. 

. Reducing the effects of beam B by coating the outer cylinder 
surface with a thin film which causes destructive interference 
of B in the coating. 

. Increasing the intensity I by coating the interior wall of 

the cylinder and increasing the intensity 1^ with a coating 
on the elastomer. 

For the present work, none of the above techniques is entirely satis 
factory. It is desirable, for example, not to have to coat either 
surface of the cylinder. Such coatings are expensive and the one on 
the inner surface is prone to scratching and wear. It is also de- 
sirable not to coat the seal surface with a material which affects 
its deformation characteristics. Finally, it is very desirable that 
the experimental apparatus employ production seals rather than seals 
especially produced with optically smooth surfaces. The cost for the 
specially produced seals can be very high, their behavior may not be 
typical of production seals, and the optically smooth surfaces can 
be prone to scratching during testing. 

Experimental work conducted during the course of the present effort 
resulted in a technique which gives good interference fringes and 
which is appropriate to the experimental apparatus. Interference 
fringes were observed and photographed with a laboratory bench setup 
by using: 

. A Parker Viton T Seal which was coated with shiny black lacquer. 
The lacquer was employed because it is flexible (to conform to 
the deformed rubber shape), strongly adherent to the rubber, 
and insoluble in the production oil used. Additionally, the 
lacquer coating does not crack easily or wear rapidly. It 
can be applied such that the thickness is small and relatively 
uniform. 

. An optical flat having a refractive index in the range of 1.5 
to 1.6. The optical flat could be moved manually to simulate 
the motion of the glass cylinder. 

. An automatic transmission oil which is a candidate for use in 
actual Stirling engines. This oil (Mobil XRL 1032 AR-1) has 
a refractive index of 1.4605. 

. A fiber optics white light source having a high intensity. 

The end of the fiber optics tube was mounted close to the op- 
tical flat and oriented so that the angle between the incident 
light A and the normal to the surface of the optical flat was 

about 30°. This mounting method eliminates the effects of 
beam B by taking advantage of the (axial) distance between 
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beams C and B; i.e., when sighting along C and D, the beam B 
does not enter the field of view. 

. A stereo microscope with a maximum magnification of about 50X. 

The fringes obtained were of sufficient optical quality that 35 mm 
color slides could be taken and a 16 mm color movie could be made. 

A black and white photograph made from a 35mm slide is shown in 
Figure 19. The dark zones are interference fringes which appear as 
multicolor bands in the original slide. The bands are very apparent 
in the color movies. The bands show the presence of a lubricant film 
between the glass and rubber at relatively low sliding speeds (less 
than 0.25 m/s (10 inches/second)) and at relatively large 
seal preloads (contact width of the seal approximately equal to half 
the seal width) . 

Once good quality interference fringes had been produced using the 
laboratory bench setup, additional optical work was postponed until 
completion of the experimental apparatus. The optical instrumentation 
was then installed on the device as shown in Figure 20. In the figure, 
the fiber optics tube can be seen at the lower center of the photo- 
graph. The microscope is at the left center. 

A low magnification photograph, taken through the microscope is 
given as Figure 21. The photograph was taken with the cylinder station- 
ary. The figure shows the aluminum bands on the test seal holder, 
the two fiber backing rings for the T seal, and the T seal itself. 

At this low magnification, the contact zone cannot be seen but is 
approximately equidistant from the two fiber backing rings. A photo- 
graph for which a higher magnification was used is given as Figure 22. 
For this figure, a scale was taped to the outside of the transparent 
cylinder. Each of the smallest divisions is one millimeter. In 
Figure 22, the contact zone can be seen bounded by horizontal light 
lines in the dark region of the T seal. 

Interferometry was attempted with both a stationary transparent 
cylinder and for a slowly moving transparent cylinder (the flywheel 
was turned by hand). Fringes were obtained, but were not sufficient 
quality to be photographed. The lack of optical quality was caused 
by the transparent Lexan cylinder, whose optical characteristics were 
poor. Specifically, both the inner and outer surfaces of the cylinder 
had circumferential scratches and grooves which appeared to have been 
produced by turning during manufacture. These scratches and grooves 
tended to dissipate both the incident and emerging light. In addition, 
the contact region, during motion of the cylinder, was rather irregular 
on a microscopic level. This caused the fringes to be numerous, 
close together, and non-uniform. 

An attempt was made to polish the Lexan cylinder. This did improve 
the surface quality somewhat, but not sufficiently that a smooth 
contact zone resulted and that sufficient optical quality was attained. 
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Figure 19 Interference Fringes Produced with Laboratory 
Bench Setup (Original Color Photograph Made 
Using White Light) 
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Figure 20 Experimental Apparatus With 

Optical Instrumentation Installed 
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Figure 22 Close-up of "T" Seal in Contact 
With Transparent Cylinder 
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The attempt at polishing did indicate that Lexan does not polish well 
because of its inherent softness. In addition, the material does not 
maintain its degree of polish well — again due to its inherent softness. 

The above work led to the conclusion that Lexan is unsuitable for use 
in the transparent cylinder. Two alternative materials are glass and 
acrylic. At the time of this writing, a glass cylinder is on order 
and due for installation in the experimental device before the end 
of March. The manufacturing of the glass, however, is risky and 
it is not certain that a glass cylinder having the proper dimensions 
and tolerances can be secured. Nevertheless, the glass cylinder is 
the most desirable, since the inner and outer surfaces will be op- 
tically smooth and will be very scratch resistant. 

Acrylic is the second choice — it is desirable because it is relatively 
easy to polish and relatively scratch resistant. In the event that 
substantial problems arise in manufacturing the glass cylinder, this 
material will be used in the apparatus. However, since acrylic is 
susceptable to scratching, designs are being considered currently in 
which the acrylic member is expendable and easily replaced. 
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ANALYTICAL REPORT 


In this section, the analytical work completed during the past year is 
discussed. The analytical work was directed at determining the film 
thickness distribution at the seal/cylinder interface. This film 
thickness distribution is important because it controls the net leakage 
during each cycle of the reciprocating motion. The film thickness 
distribution also affects the frictional forces produced by the 
sliding seal. 

Determination of the film thickness distribution for a reciprocating 
elastomeric seal requires that two elements of seal behavior be 
considered. Firstly, due to the low modulus of the seal, the fllm 
thickness and the pressure distributions are coupled. These ms be 
determined together by theories developed in elastohydrodynamic (EH ) 
lubrication. Secondly, the film thickness at the seal/cylmder inter- 
face fluctuates during each reciprocating cycle. This is unlike con 
ventional EHD problems for bearings and cams where steady state EHD 
theories provide accurate film thickness predictions. Because of the 
reciprocating motion, the film thickness is influenced by both squeeze 
film and sliding effects - an accurate determination of film thickness 
behavior therefore requires a solution of the time-dependent Reynolds 
equation coupled with the elasticity equation. However, a reasonable 
estimate of the average film thickness level in rod seals can be had. 
by using the existing steady-state EHD theories based on some effective 
velocity over each cycle of the reciprocating motion. Contributions 
in this category include the earlier analysis by Dowson and Swales 
[4, 5], and later by Hooke, et. al. [6, 7]. 


The analytical work given in this section includes two parts. First, 
an approximate method is presented to estimate the mid-stroke and the 
end stroke film thickness. In the second part, a time-dependent EHD 
analysis is developed to determine the temporal behavior of the film 
thickness distribution during the entire cycle of the reciprocating 
motion. The nomenclature used for both parts of the analytical wor 
is given in Table 1. 


Es t imate of the Mid-Stroke and End-Stroke Film Thickness 


The variation of film profile in a rod seal during a reciprocating cycle 
was observed first by Blok and Koens [8] using the interferometric 
technique. Figure 23 shows the observed film profiles at twelve 

_7V 

different stroke positions. At mid-stroke positions, i.e., ut - ^ * 


ut = the sliding velocity is at its maximum and the film profile 

is similar to those predicted by the steady-state EHD analyses. At 

^ 1 _ j_1 — 1 , 


the end-stroke positions, i.e. 


tot = 0 and wt = it, the sliding velocity 


is zero, and surfaces are separated only by squeeze-film action. As 
the shaft travels from mid-stroke to end-stroke position, the minimum 
film shifts forward from the exit side of the contact to the entrance 
region. There also seems to be a considerable variation of the average 
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Table 1 

Nomenclature for Analytical Report 

a = semi-major axis of the contact ellipse 

b = semi-minor axis of the contact ellipse 

c = an integration constant used in Equation (32) 

c' = an integration constant used in Equation (36) 

c ' , c ' = integration constant c' in inlet and exit half of the 

° contact region respectively 



d^ dg k = derivative coefficients for uneven grids 

E^, E£ = elastic moduli for bodies 1 and 2 



f = frequency of reciprocating cycle 

h = film thickness 

h = center film thickness 
o 

h , = minimum film thickness 

min 


h 


o 


h 

o 

R 


h T = h for the last time interval or last stroke position 
o,L o 


H 


h 

h 


H* 


H 

k,L 


1/2 

, used by Herrebrugh [9] 

H at kth grid 

H, for the last stroke position 

K 


“min E'r I 

R 2 p ul 

o 1 
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j = grid designation 

k = grid designation, also used as the ellipticity ratio by 
Hamrock and Dowson [10] 

k Q = grid number at the contact center 
k ff = grid number at the film termination 


M = 


(e^r) » used by Herrebru S h t 9 ] 


p = pressure 

p = pressure at the contact center or maximum Hertzian pressure 

P = *- 
P o 

P = P at kth grid 

K 

0, . = a kernel function, see Appendix B 

k,j 

R = radius of "0" ring 

R , R = effective radius in the x and y direction respectively 
x y 

S = stroke of the shaft 

t = time 


T = 


u t 
o 


u = entrainment velocity 


0)S 


u q = entrainment velocity at mid-stroke - — 


U = 


y u 
o o 

E'R 


V = normal approach velocity 

W = load per unit circumferential length of "0" ring 
W^, = total load 


W„ 


W = 


E'R 

x 


x = coordinate in the direction of motion 
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cr|x 


inlet coordinate 


= coordinate at the film termination 
viscosity 

= Poisson's ratio of bodies 1 and 2 
2nf 

equations to be solved for P^» etc., see Equations (37) 
dummy coordinate for x 



film thickness during the reciprocating cycle. It will be seen later 
in the detailed time-dependent EHD analysis that these observed film 
characteristics are completely in accord with the predicted film pro- 
files. It should be noted that if the kinematics are the same during 
the forward and backward stroke, the film profiles should also have 
the same characteristics. Figure 24 depicts the qualitative features 
of the film profiles relative to the four stroke positions in Stirling 
cycle rod seals. 

Approximate Mid-Stroke Film Predictions. Even though results of steady- 
state EHD analyses are not exactly applicable to a reciprocating rod 
seal, they are nevertheless suitable for making approximate predictions 
of the mid-stroke film thickness and the level of film thickness for 
each cycle. For this purpose, results obtained by Herrebrugh [ 9 ] for 
line contacts as well as those recently developed by Hamrock and Dowson 
[ 10 ] for low modulus point contacts have been adopted. 

The numerical solution by Herrebrugh, Figure 25, gives the variation 
- of a dimensionless film parameter H* with an EHD parameter M. Since 
the load parameter, W/ER, in rod seals is usually very large, the 
portion of the curve in Figure 25 applicable to rod seals is towards 
right in the heavily loaded region. In this region, the curve can be 
best fitted by 

H* = 2.45 M _1/ ^ (1) 

where 



Equation (1) can be rearranged to give 



Since in line contacts 



(4) 


( 5 ) 
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Figure 25 The Functional Relationship: H* = F(M). 

Dash and dot Gilmbe 1-Mar tin; 
Dashed line — Weber and Saalfeld; 
Circled points — author's results 


Equation (4) in terms of a pressure parameter -gr , becomes 



An alternate method can also be used to estimate the minimum film at 
the mid-stroke by modifying the minimum film thickness formula given 
recently by Hamrock and Dowson [10] for low-modulus point contacts. 
Reference [10] gives 


- 7.43 (1 - 0.85 e-°- 31k ) U 0 ' 65 w'°' 2 ‘ 
R 

x 


where 


h . = 

mm 


a 

b 

R 

ji 

R 

y 

u 

v 


minimum film thickness 
0.64 


" t ' 1 - 03 (i*) 

» V * 


X* 

semi-major axis of contact ellipse 
semi-minor axis of contact ellipse 
effective radius in the direction of rolling 
effective radius perpendicular to rolling 


E'R 


viscosity 

entrainment velocity 


1 

E' 


W 


+ t^l) 

2 l \ E 2 j 


W„ 


E'R 


X 


(7) 
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W T = total load of the contact 


From Hertz theory for an elliptical contact, it can be shown that 


o _ _£ 

E' 2 it 


W_ 


x * 


2 ' l E'R 2 I 

X 


( 8 ) 


or 


W„ 


E’R 2 
x 


2i[ a b 
3 b R 


-l 2 ^ 

J E 


(9) 


Substitution of [9] into [7] gives 


h . 
nun 

R 

x 


= 7.43 (1 - 0.85 


-0.31k. 0.65 

e ) U 


,,v-.2lf2ir P o fb \ 2 1 

(k) [-3 r (irJ J 


- 0.21 


( 10 ) 


It was suggested in Reference [10] that for low-modulus contact 
Equation (10) can be converted to line contact if k is taken to be 11. 
Setting k = 11 and using the following relation for line contacts. 


b^ 

R 



( 11 ) 


Equation (10), for line contacts, becomes 


h - • 63 

m±n = 7.43 (1 - 0.85 e -0,31 x U ) (ll)“* 21 (16--^) _ * 21 u°* 65 (-^2_) 


R 


x 


( 12 ) 


or 


- - mj - n = 2 09 — 

R U E’ 

x 


,-.63 


(13) 


Equation (13) is the minimum film for low-modulus line contacts based 
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on Hamrock and Dowson's point contact theory. Film thickness curves 
based on Equation (6) and Equation (13) for conditions used in the 
Stirling engine rod seals are plotted in Figure 26. The agreement 
between these two theories is seen to be reasonably good at low loads. 
At higher loads, the larger exponent in Equation (13), deduced from 
the point contact theory, yields a much lower estimate compared to 
Herrebrugh's theory. Nevertheless, for pressures encountered in 

Stirling engine seals 3.5 x 1*"* • 14 x 10 kPa (500—2000 psi) 

the discrepancy is not serious; therefore, both theories appear to be 

satisfactory . 

Herrebrugh's theory was recently correlated experimentaly by Swales 
et. al. with rubber pads [5]. Comparison in Reference [5] showed that 
the load dependence in the experimental data is considerably smaller 
than in Herrebrugh's formula. Since the values of P q /E i- n Swales 

experiment were considerably smaller than would be encountered in 
Stirling engine seals, it is believed that film predictions based 
on Equation (6) or Equation (13) for rod seals are likely to be on 
the pessimistic side. 

End-Stroke Film Thickness . Based on the steady state analysis, the 
film thickness at both ends of the stroke would be zero because the 
entrainment velocity is zero. However, due to squeeze-film effects, 
the end-stroke film thickness will not vanish, but will attain a finite 
value which depends on the squeeze-film velocity. 

To estimate the values of end-stroke film thickness due to squeeze- 
film effects, the results developed by Herrebrugh [ll] for two normally 
approaching cylinders with an isoviscous lubricant are used. Based 
on his findings, the center film thickness during a normal approach of 
two cylinders in the heavily loaded region is 



where 


V = normal approaching velocity 

It is assumed that the film thickness at the end-stroke does not deviate 
greatly from the mid-stroke film thickness; i.e.. 


"o 


(15) 
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(microinch) 



* H (p*D 


Figure 26 Comparison of Approximate Analytical 
Film Thickness Results 
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The normal approaching velocity can then be approximated by 

h - h 

Ah o °’f ° >7r 

V = = to (16) 

where the time increment can be taken as one quarter of the period; i.e., 

* - i ID < 17) 

In Equation (17), f is the frequency of the strokes. 


Rearranging of Equation (14) gives 

3 


u V |h V 
H o _ I o 1 

E'R \R / 


1 


(4.85)- 


Ah, 




(18) 


By substituting — — for V one obtains 



(19) 


Here, u is the entrainment velocity at the mid-stroke, 
o 


u 

o 



1 

2 


from which 


irfS/2 


u At 
o 



1 

4f 


ttS 

8 


so that 


( 20 ) 
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Equation ( 6 ) can now be used to express h c /R in terms of the speed and 
load parameter. Substituting Equations (20) and ( 6 ) into (19) yields 


~TT = (if) 2 * 385 x u °’ 25 v a ' c 1 g (-t) 

h V 8 R l 285.8 l E' I 

o * 9 

=0.00327 (|) x U 0 - 25 (^J 

S P Q 

When the following ranges for U, and —7 are used, 

K E 

! - 4 - 12 


( 21 ) 


( 22 ) 


U = 2 x 10 6 - 5 x 10 7 
p 

ttt = 0.05 - 0.2 

E 

the ratio ~r~ is found to be within the range from 0.59 - 0.0087 
h 
o 

This shows that the squeeze-film effect can be quite substantial under 

conditions in rod seal lubrication. The variation of the center film 

h during the entire cycle of the reciprocating motion is expected to 
c 

be small for high values of p /e'. 

Transient Film Thickness Analysis 

It was shown in the last section that the mid-stroke film thickness can 
be estimated by equations derived from existing low-modulus EHD theories. 
The reduction of film thickness at the end-stroke can be estimated by 
a relation deduced from Herrebrugh's normal approach EHD analysis. In 
this section, an analysis is developed to determine the temporal behavior 
of the film thickness by solving simultaneously the time-dependent 
Reynolds' equation and the deformation equation. Since it was antici- 
pated that this transient EHD solution for severly loaded rod seal 
contacts might encounter unforseen difficulties, it was decided that 
the first attempt at producing the transient analysis be based on simple 
Hertzian contact. The analysis below, therefore, does not include 
the deformation due to contact pressures exerted by groove surfaces on 
the seal ring. The analytical effort was focused on developing a 
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numerical method to handle the transient EHD solution under severe 
p^gggiir’Q loading# Other deformation profiles, including those produced 
by groove contact pressures, can be incorporated into the perfected 
method at a later date. 


Governing Equations. The time-dependent Reynolds equation for a pure 
Hertzian contact is 


_ 9 _ 

9x 


(» 3 l) ' 12 "o ( 


, . . 9h . 9h 
= 12y o I u(t) ^ + "9tJ 


(23) 


The film thickness is governed by 


h = h +^r 
o 2 R 


«/ X. 


ff 


p(C) In 


g -x 


dg 


(24) 


where 

x - coordinate along the sliding surface 

h = film thickness 

p = pressure 

u = viscosity 

o 

u = time varying velocity 

h = center film at x = 0 

o 

R = seal radius 

x. = inlet boundary of the lubricant film 
= exit boundary of the lubricant film 
5 = dummy variable for x 


Introducing the following non-dimensional quantities: 

- x 

X = b 

b = half contact width 


o 

P = p/p Q 

p = maximum Hertzian pressure 

r rt 
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h 

o 


(v r ) 


V 


48 Vo 
E'R 


U o “ mid-stroke entrainment velocity, half of the mid-stroke sliding 
velocity ( u = wS/4) 

C = 16 (p /E') 2 /(h /R) 

± o o 

S = stroke of the piston 

rp _ u t 

T - o 

b 

to = angular velocity = 27rf 

l = C/b 


Equation (23) and Equation (24) in non-dimensional form become 


3 _ 

3x 




s intot 



3(h H)1 
o 

3T 


H = 1 + C 


-2 

x 

2 



P(£) In dt 

\l\ 


(25) 


(26) 


In the above, the shaft velocity is assumed to be sinusoidal. The 
boundary conditions for Equation (25) are: 


P = 0 at x = x. 

x 


P = 0 and — =0 at x = x 

dx “ 


P = 1 at x = 0 



(27) 

(28) 

(29) 

(30) 


For a given set of parameters ° , and , one can determine the 

El El K K 

center film thickness and the film thickness and pressure profiles 
at each time interval using numerical schemes. The two numerical 
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schemes considered in solving Equation (25) and (26) and the experience 
gathered with each are documented below. 

Numerical Treatment — Direct Discretization (Scheme "A") . This scheme 
is the conventional way of handling the Reynolds' equation. In 
Equation (25), the second order derivative of P is directly discretized 
at a grid point k. The discretized equation appears as 


\ = 


- H. 


k-1/2 Ax 


^k ^k-1 3 

k — + H 


P -P 
k+1 k 


k-1 


- smrnt 


k+1/ 2 Ax, 


h 2 
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[ d 3,k \+l + d 2,k H k d l,k H k-lj 
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where 


H 
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1 + C, 




P. 
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C x = 16( P o/E') 2 / (^) 
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(31) 


P. 

3 


= x coordinate at grid 
= P at grid j 


k 


Q 


kj 


= A Kernel 
duced in 


function developed in Reference [ 10] , and repro- 
Appendix B 
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l,k 




Ax k-l (Ax k + A^) 


2,k 


A \ ~ A Vl 
A *k Ax k-1 


3,k 


Ax. 


k-1 


4; k ( 4 ; k + 4^) 


h T = h for the previous 
o, L o 


H, = H, for the previous 

K y J_» K 


AT 


u At 
o 


time interval 


time interval 


At each time step, Equation (31) is a set of equations written at grids 
2 to k ff -l for a total of k ff -2 equations. These k ff -2 equations can 

be solved by the Newton-Raphson procedure to obtain values for k^-2 


unknowns consisting of to ? k ^* k to ^ k an d h Q . 

. O o ' L ff 

solution. 


In this 


P 


1 


0 


P k =1 
o 


P 


v ff 


0 


which are the boundary conditions for Equation (31) . 

The location of x, is not fixed, and must be such to make the slope 
dP _ ff 

_ vanish at x, . This is accomplished by moving x, towards the 
dx k ff k ff 

contact center one grid point at each iteration whenever a negative 
pressure appears near k^. 
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Scheme "A" was tried and failed to yield satisfactory results particu- 
larly at high loads. It begins to experience difficulties at 

Z 2 . = 0.01. It is not known whether the difficulties arose because of 

the inherent discretization error associated with the interpolations 

factors d. . , d, , , d for uneven grids or because of a possible 

ljlc Zjk. JjK- t , , 

error in the complicated algebraic expression of k which Has to be 

3h 

o 

evaluated in the Newton-Raphson procedure. Scheme " A " M was fi n all y 
abandoned in favor of an alternative approach. Scheme "B' , which is 
described next. 

Numeric al Treatment — Discretization of Integrated Reyno lds Equatio n. 
(Scheme "B") . This scheme follows closely that developed previously 
by Lee and Cheng [12] for studying the effect of a single asperity 
entering the inlet zone of an EHD contact. In this formulation, the 
Reynolds equation is first integrated once with respect to R, and then 

discretized. 


Integrating Equation (25) once, one obtains 



H~ 


3x 


smut 


"*irf 

O X. 

1 


x 3(h H) 
o 

3T 


d£ + c 


(32) 


where c is an integration constant. 


If is zero at x = 0, then 
dx 


c = - sinmt - 


tf 

X . 

1 


8(hH) 
3T 


dt 


(33) 


For general cases, at x = 0 is. small but not exactly zero; therefore, 

a small residual constant c ' must be provided to permit any small non- 
zero pressure gradient to exist at x = 0. Thus, for general cases. 


c = c ' - sitiyt - — 


i-r 

° ^x. 

l 


8 de 

3T 


(34) 


Substituting (34) in (32), and observing that 
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(35) 
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Equation (32) becomes 


h 2 H 3 — = sincjt (H - 1) - — 

° 

8X iln W 
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1. f° 3 < 5 o H > 

hJ- 8T 
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(36) 


Equation (36) is discretized at the mid-point, k + 1/2, between two 
pressure nodes, and The discretized equation is 


* , R 2 [ 3 Vl ~ P k 
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By referring to the following diagram, 



Equations (37) can be written at k = 1 to k = k^ - 1 for a total 
number of k^ - 1 equations. In solving Equations (37), the following 
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pressures are known, 



Thus, the total number of unknown pressure nodes is k^ f - 3, and they 

are P„ to P, , and P, , , to P. , . These k,. -3 unknown pressure nodes 
2 k -1 k +1 k^-1 it 

o o f f 

together with h and c' form the right number of unknowns for the 
k^-1 equations. 


As mentioned in Scheme "A", must satisfy the free boundary condition 

dP 

— =0. This is accomplished at each iteration by searching for the 

dx 

point near the exit region where the following term vanishes 


sinrnt (H - 1) - — 


i r° 9( v 

h j- 9t 

n V 


It was found the inclusion of h as an unknown in Equations (37) in 

o 

the Newton-Raphson procedure gave rise to some difficulties in conver- 
gence. This was remedied by allowing the integration constant c' in 
the inlet half to be different from the value of c' in the exit half. 

In so doing, c' is now replaced with two unknown constants, c.j. for the 

inlet region, and c' for the exit region; however, h Q is now consid- 
ered as known. Thus, the total number of unknowns still matches the 
total number of equations. 

At each time step, one begins the first iteration with an initial guess 

of h , which is usually the value of h for the previous time step, 
o _ o 

With the first guess of h Q , Equations (37) can be solved for P^ and C.J. 

and c ' If c' is not equal to c', h is readjusted until c' matches c'. 
o*I _ oo I o 

This procedure permits h Q values to be calculated during half of the 
reciprocating cycle. However, h Q thus determined in the first half 
cycle may not satisfy the periodic conditions, i.e., 
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^ h o^u)t=iT/2 wt=3ir/2 

The cycle must be repeated until the periodic boundary condition in time 
is satisfied. 

The procedure seems to work well for a wide range of load and speed 
parameters, even for conditions beyond those occuring in Stirling engine 
rod seals. Results obtained from use of the procedure are given in the 
next section. In the Appendicies, details and use of the procedure^ 
are discussed. Specifically, the expressions for the derivatives of 

ill with respect to P. for solving Equations (37) can be found in 

/ 2 j 

Appendix C. The computer program for Scheme "B" is entitled "RODSLJ" 
and is documented in Appendix A. 


Results 

Three input parameters are required for each run. 
P 

— r : load parameter 


These include: 


p u 

U = : speed parameter 

E R 


— : stroke to seal radius ratio 

R 


For the conditions to be used in the experimental rig, the range of the 
above dimensionless parameters was found to be: 


p 

— = 0.05 to 0.25 

U = 10 -7 to 2.0 x 10' 


| = 4 to 16 

Table 2 lists the input data for seven runs. Runs 50, 51, and 52 are 
for finding the effect of speed; Runs 51, 53, and 54 for the effect of 
load; and Runs 51, 55, and 56 for the effect of stroke to radius ratio. 
Figures 27 to 34 show dimensionless film profiles at six different 
stroke positions during the forward stroke (u)t=0 to cot=TT) of the re- 
ciprocating cycle. Since the sealing pressure differential is not 
considered in the present analysis, the film profiles during the back- 
ward stroke (uit=ir to u>t=2Tr) should be a mirror image of the forward 
stroke profiles, which are the only ones presented here. All film 
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Table 2 List of Input Data 


Run No. 

P o 

E' 

U 

S 

R 

y-v- 

50 

0.1 

1 

o 

tH 

8 

31.6 

51 

0.1 

5xl0“ 7 

8 

14.14 

52 

0.1 

2x10" 6 

8 

7.07 

53 

0.25 

5xl0~ ? 

8 

88.3 

54 

0.05 

5x10" 7 

8 

3.54 

55 

0.1 

5x1 0" 7 

16 

14.14 

56 

0.1 

5x10“ 7 

4 

14.14 
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profiles are normalized with respect to the film thickness at the center. 

The value of h , the non-dimensional film thickness at the center, at 
o 

each stroke position is labelled for each profile. 

For Run 54, the most lightly loaded case, the film profile between the 

last half of the forward stroke (wt = ^ to ut =ir) and the first half of 

the backward stroke (wt =tt to wt = are replotted on the same scale 

to show their relative positions. These curves show a trend consistent 
with these observed by Blok and Koens [8]. It would be of interest to 
ascertain the exact conditions used in their experimental runs so that 
more quantitative comparisons can be made. 

Aside from the fact that the edge protrusions are much more pronounced 
for heavily loaded or low speed cases, the general features of film 
profiles are very similar for all runs. It is interesting to note that 
at the end-stroke position (mt = tt), where the sliding velocity vanishes, 
the film droops down more in the inlet half region, whereas in the exit 
half the film remains reasonably flat or even rises slightly above the 
center film in some cases. This feature seems to be due to the sliding 
action producing a very effective suction once a negative wedge is 
formed at the inlet. As can be seen in Figure 23, this same feature 
also exists in the observed film profiles at the end-strokes. 

The variation of the center film thickness during half of the recipro- 
cating cycle (wt = ~ bo tot = ^ ) is shown in Figure 35. At the mid- 

stroke (wt = ^), the shaft velocity reaches a maximum and is decreasing. 

However, the center film thickness does not reach a maximum until about 
10-15 degrees after the mid-stroke. This shift is due to the fact that 
part of the sliding action at the mid-stroke is nullified by the 

squeeze— film action "“r. The stronger the squeeze— film term, the larger 

0 J. 

the shift of the maximum film from the mid-stroke position. As the 
shaft reaches the end-stroke position (wt =ir) , the sliding action ceases, 
and the pressure within the film is supported entirely by the squeeze- 
film action. Perhaps the most interesting feature in this curve is the 

sudden reduction of h soon after the end— stroke position, where the 

o 

motion reversal takes place. This reduction is attributable to the pen- 
etration of the exit protrusion into the center of the contact after 
the shaft reverses its direction. This penetration causes a sharp re- 
duction h for low-speed or high-pressure cases. It becomes more pro- 
o 

nounced as the squeeze-film effect becomes more influential. 

Figure 36 shows influence of speed on h Q . It is seen that the squeeze- 
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0)t 


Figure 35 Variation of h with the Stroke Position 

o 


Run #50, 51 and 52 




film a; tion is not altered significantly even for a five-fold increase 
in speed. This trend agrees well with the approximate analysis of 
squeeze-film effect represented by Equation (20) . 


Figure 37 shows the effect of load parameter on h Q . It is seen that 

the fluctuation of h is most strongly influenced by the change in 

o 

pressure. A five-fold increase in Hertzian pressure from 0.05 to 0.25 


yields more than ten times reduction of h Q during the last half of 
the forward cycle (ait = to ait =ir) . This result again confirms 


qualitatively the approximate formula Equation (20) derived from 
Herrebrugh's squeeze-film analysis. It should be noted that for the 


highest Hertzian pressure 



.25 


the reduction of h 

o 


right after 


the end-stroke position (ait =tt) is quite sudden. It should also be 
noted that the lowest film thickness at this position i s accompanied 
by a substantial sliding velocity. If the seal surface is to fail by 
sliding wear, it is likely to begin first at this position. 


Figure 38 shows the effect of stroke length on the variations of h Q . 

As expected 5 a longer stroke would lessen the squeeze-film effect 

and allow h to approach more to the steady-state value based solely 
o 

on the sliding velocity. With a shorter stroke, the squeeze-film 
action dominates, and the fluctuation due to the sinusoidal sliding 
velocity is considerably suppressed. 


Concluding Remarks on Analyses 

A time-dependent numerical EHD analysis for low-modulus line contacts 
was developed to investigate the temporal behavior of the film pro- 
files of sliding rod seals during a reciprocating cycle. Resulting 
film profiles reveal all essential features of the temporal film 
behavior observed by Blok and Koens [8] between a rubber cylinder and 
a transparent reciprocating surface. 

The results indicate that the squeeze-film effect dampens the film 
thickness fluctuations caused by the oscillatory sliding velocity. It 
also produces a shift of the maximum film from the mid-stroke position, 
where the maximum velocity occurs. Likewise, the minimum center film 
thickness occurs not at the end-stroke, but a few degrees thereafter. 
For cases of pronounced squeeze-film effect, usually associated with 
high-load cases, the reduction of the center film after the end- 
stroke is rather abrupt, particularly at a position where there is 
substantial sliding velocity. 
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The squeeze-film effect increases strongly with an increase in contact 
pressure, but is almost uninfluenced by a significant change in speed. 
These findings correlate well with an approximate formula for the 
change of film thickness due to squeeze-film action based on an existing 
normal approach analysis [ll]. 

The results developed in the present analysis are based on a pure 
Hertzian contact uninfluenced by any contact pressure between the seal 
ring and the groove walls or by the effects of a large sealing pressure 
differential. The analysis should be valid for cases where the sealing 
pressure differential is small compared to the contact pressure. The 
transient analysis developed is an essential building block for further 
study of the temporal film behavior of rod seals under large fluid 
sealing pressures. The analysis can be readily modified to include an 
inlet deformation profile obtained from a static finite element defor- 
mation computation. When such a profile is combined with the present 
analysis, the resulting analysis will be capable of treating the wall 
pressure effects and the pressure drop across the contact region. 
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CONCLUSIONS AND RECOMMENDATIONS 


On the basis of work done to date, the following conclusions have 
been reached regarding the experimental apparatus and the analytical 
effort: 

Experimental Apparatus 

. The experimental apparatus, in design and performance, is 
suitable for evaluating Stirling engine rod seal behavior . 

The apparatus produces a well-defined reciprocating motion 
of the transparent cylinder. It exhibits little vibration. 
The built-in force cell measures properly axial forces 
applied to the test seal holder. In addition, the device 
is rugged and easy to use. 

The interferometric measurements, necessary in the evaluation 
of seal behavior, can be made with the device. However, the 
Lexan presently used for the transparent cylinder must be re- 
placed. The new material should be optically smooth both on 
the inner and outer diameters and should be scratch resis- 
tant. The first choice for this material is glass and the 
second choice is acrylic. 


Analysis 

. A time dependent numerical method was developed which deter- 
mines the fluid film thickness and pressure distribution at 
the contact zone between an elastomeric seal and a recipro- 
cating cylinder. The method is appropriate for the high 
values of static contact pressure to seal modulus ratio found 
in Stirling engine rod seal applications. At its present 
stage of development, the method applies for a Hertzian static 
contact pressure profile and for zero pressure gradient across 
the seal. 

. Application of the numerical method produces results which 
correlate well with experimental results in the literature 
and with appropriate approximate analyses. The results 
obtained show that 

a) the squeeze-film effect dampens the film thickness fluc- 
tuations caused by the oscillatory sliding velocity 

b) the maximum film does not occur at the instant of max- 
imum velocity 

c) the minimum film does not occur at the instant of zero 
velocity 

d) the squeeze-film effect is affected strongly by the 
contact pressure but is almost independent of the sliding 
velocity. 
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Continuation of the contract work should be directed at maximizing 
the usefulness of the analytical and experimental tools developed to 
date. The following recommendations are based on this objective. 

It is recommended that the experimental apparatus be used in an in- 
depth study of Stirling engine rod seal behavior. To do this requires 
that some refinement to the apparatus be done, that a test program 
be developed and conducted, that the resulting data be reduced ap- 
propriately, and that some runs of the analytical seal model be per- 
formed for cases corresponding to the experimental situations. 

The refinement of the experimental apparatus entails the addition of 
hardware to synchronize a stroboscopic photographic flash with the 
cylinder position. 

The test program should include use of the apparatus with the trans- 
parent cylinder, with an aluminum cylinder, and with a coated and 
uncoated seal. Using the transparent cylinder and the coated seal 
can provide, primarily, film thickness data for comparison with theo- 
retical results. Parameters which can be considered include: 

1. Seal type and seal material 

2. Seal preload and seal groove geometry 

3. Gas pressure to 690 kPa (100 psi), (zero pressure 
for at least one set of tests) 

4. Frequency (10 to 50 Hz) 

5. Temperature and lubricant type 

6. Seal surface finish (variations to be obtained by changing 
the optical coating on the seal) 

Operation of the apparatus with the glass cylinder and an uncoated 
seal can reveal the role of the optical coating in the results with 
the coated seal. To accomplish this, measurements of leakage and 
friction-force can be made and compared to these obtained for the 
coated seal. 

Operation of the apparatus with the aluminum cylinder and an uncoated 
seal can reveal the role of cylinder surface finish on friction and 
leakage. In this, measurements of friction and leakage can be taken 
and compared to those obtained for the uncoated seal/glass cylinder 
combination. Also, aluminum cylinders having various surface finishes 
can be employed in the comparisons. 

Regarding data reduction, the most important aspect is quantifying 
the interferometry results to determine the film thickness and the 
film thickness profile. Techniques for the film thickness measure- 
ment include white light with color tables, two discrete wavelengths 
of light, or monochromatic light with counting of fringes born. It 
is desirable that the most suitable technique or techniques be deter- 
mined and used in the above experimental work. 
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Additional runs of the elastohydrodynamic seal model will be needed 
for those situations to be treated experimentally. Such situations, 
not yet considered analytically, include various combinations of seal 
properties, seal preload, cyclic frequency, and lubricant viscosity. 
Comparing similar experinental/analytical cases (to the extent that 
the model can represent the experimental situation) will provide 
insight into both the important aspects of the elastohydrodynamic seal 
behavior and additional capabilities desirable for eventual addition 
to the computer model. Capabilities already identified for eventual 
addition to the computer model are a general (non-Hertzian) static 
pressure distribution and a pressure drop across the seal. 
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APPENDIX A 


COMPUTER CODE "RODSLJ" 



APPENDIX A 


Computer Code "RQDSLJ" 


Flow Chart 

The computer code "RODSLJ" was written to calculate the time-dependent 
pressure and film thickness profiles, P^ and H^, at various stroke 

positions during half of the reciprocating cycle (mt=7r/2 to o)t=3ir/2) . 

The flow diagram shown in Figure A— 1 is based on the numerical treatment 
described as Scheme "B" in the Analytical Report. For a set of given 
input data, the calculation of P fc and for the entire half cycle is 

nested in four loops. In the innermost loop (IT=1, ITP) , P^ and are 

calculated for a prescribed h Q to satisfy Equations (37) . The second 

inner loop (ITH=1, ITH) is an iterative loop for solving the center 

film h to satisfy the condition, c ' - c * =0. The third inner loop 

(MT=MII, MTMAX) performs the calculation of h Q , P k> and H k in equal 

time intervals from the one mid-stroke position ( tot=Tr / 2 ) to the other 
mid-stroke position (o)t=3ir/2) . The outermost loop (IE=I, ITE) repeats 

the MT loop until h ,„ = h „ 

o,tt/2 o,37t/2 

The nomenclature used in the computer program is listed in Table A-l. 














Figure A-l Flow Diagram for Code "RODSLJ" (Continued from previous page) 





Table A-l 

List of Major Computer Symbols 
Computer Symbol Symbol in the Analysis 


KA 

k A 

K4> 

k 

o 

KF 

k f 

KFF 

k ff 

DX(K) 

A *k 

X(K) 

*k 

P <f)E 

p /E' 
o 

UB 

U 

WT 

0)t 

DWT 

Acot 

FL4B 

S/ (4b) 

P(K) 

P k 

C7 

C 7 

DT 

AT 

Q(I,J) 

Q i,j 

H«frB 

h 

o 

H(K) 

\ 

CKI 

4 

CKtf) 

c’ 

o 

HL(K) 

H k,L 

H(()BL 

h _ 
o,L 

SUMS(K,TH,0) 

See Eqs. (38) 


MT 


Time interval 
designation 


Fortran Listings of Code '^ODSLJ 11 


10 

. 1.1 

20 

30 

31 

32 
40 
50 
60 
70 
80 
90 
100 
.1.03 C 

104 C 

105 C 

106 C 
1.10 C 
.1.20 C 
130 C 
140 C 
.141 
142 
150 
160 
16.1 
170 
180 
181 
190 
200 
210 
211 
220 
221 

230 

231 

240 

241 
250 
260 
270 
280 

290 

291 
300 
3 1 0 
320 
330 


P R 0 G R A M R 0 D S I... J ( I N P U T * CONS 0 L 1 0 U T P U T y T A P E 5 == INPUT y I A P E 6 -CONSOL y 
1 TAPE 7- OUT PUT ) 

DIMENSION A ( 60 y 60 ) yC(60) yWQRK(60) yl-IINT (80) ySlJMTK(SO) 

D I MENS I ON HM ( 80 ) y HP < SO ) y HM2 ( 80 ) r I-IP2 ( 80 ) y I-IM3 ( 80 ) y HP3 ( 80 ) 
DIMENSION PTEM ( 80 ) y HTEM ( 80 ) 

D I M E N S 1 0 N P S T 0 R < 2 5 y 8 0 ) y H S T 0 R ( 2 5 y 8 0 > y H 0 B S T ( 2 5 > y S U M K F F (25) 
COMMON H 0 D y H 0 B A ( 2 5 ) y W T y D T y K A y K 0 y K F y K R y K F F y li B y P 0 E 
COMMON X ( 80 ) y DX < 80 ) y Q < 80 y 80 ) y Q JA ( 80 ) 

COMMON H ( 80 > » HD < 80 ) y HS ( 80 ) y HI... ( 80 ) y HSL ( 80 ) 

COMMON P ( 80 ) y PD ( 80 ) y PS < 80 ) y DP < 80 ) y DPS ( 80 ) y PSAVE < 80 ) 

COMMON FR ( 80 ) y FS ( 80 ) y D 1 ( 80 ) y D2 ( 80 ) y D3 ( 80 ) 

COMMON TA ( 80 ) y TB ( 80 ) y TC ( 80 ) y TD ( 80 > y TE ( 80 ) y TF ( 80 ) y TO < 80 ) y TH ( 80 ) 
COMMON 00(80) y HZ (80) 

L A S T S T A T E M E N T N 0 . U S E D y 600 
LAST KFN USEDy KF5 


READ BASIC INPUT DATA 

RODSLC WITHOUT ITERATION OF HOB 


ITMP=50 

MTF =50 


NR=5 
NW=7 
NWOIJT- 7 
READ ( NR y 1 ) 

W R I T E ( N W 0 IJ T y .1. ) 

WRITE ( NWOIJT y 1 ) 

READ ( NR y 2 ) NRLJN 

DO 2000 NNR--1 y NRLJN 

R E A D ( N R y 2 ) K A y K 0 y K F y K R y K F F y N H E R R 

W R I T E ( N W 0 IJ T y 2 ) K A y KO y KF y KRy KFF 

READ ( NR y 2 ) NS 1 y NS2 y NS3 y NS4 y NS5 y NS6 y NS7 y NS8 

W R I T E < N W 0 IJ T y 2 ) N S .1. y N S 2 y N S 3 y N S 4 y N S 5 y N S 6 y N S 7 y N S 8 

READ ( NR y 2 ) I TH y I TP y I TE y MTMAX y MTMID 

WR I TE ( NWOIJT y 2 ) I TH ? I TP y I TE y MTMAX y MTM I D 

READ ( NR y 3 ) EPSH t EPSP y EPSE y X ( 1 ) 

W R I T E ( N W 0 IJ T y 1 1 ) E P 8 I I y E P S P r E P S E 
WR I TE ( NWOIJT y 1. 1 > EPSH y X ( 1 ) 

KKF-KF-1 

R E A D ( N R y 3 ) ( D X ( K ) y K 1 r K K F ) 

DO 100 K=2yKF 
1 0 0 X ( K ) = X ( K • 1 ) + D X ( K - 1 ) 

WR I TE ( NWOIJT y 5 ) ( K y X ( K ) y K= 1 r KF ) 

PI =3. 141593 

READ ( NR y 3 ) POE y UB y DWT y FL4B 
READ ( NR y 3 ) ( P ( K ) y K= 1 y KF ) 

WR I TE ( NWOIJT y 1 1 ) POE y UB y D WT y FL4B 


3 4 0 WRIT E ( N W 0 U T y 5 M K y P ( K ) y K = 1 y K F ) 

350 C7=16 . 0*P0E#*2/SQRT ( 48 . 0*UB ) 

355 DT=DWT#FL4B 


CALCULATE PD ( K ) 


360 C 

370 C 

380 C 

390 

400 

410 

420 

430 105 

440 
450 C 
460 C 
470 C 
480 
490 
500 
510 

520 106 

530 104 

540 

550 

560 

570 

580 115 

581 

583 510 

585 

586 

587 

588 

589 

590 

591 610 
630 C 
640 C 
650 C 
660 
661 
662 

663 

664 

665 

666 
667 
670 C 
680 C 
690 C 
700 
705 

710 

711 

712 

713 

714 


DO 105 K=1 y KF r 
PD<K)=0 
AXK-ABS ( X ( K ) ) 

IF ( AXK ♦ L T ♦ 1,0) PD(K)=Sa'RT(1.0-X(K>*#2) 
PS ( K ) “P ( K ) -PD ( K ) 

WR ITE ( N WOUT f 5 ) < K f PS ( K > f K= 1 1 KF ) 

CALCULATE Q(IfJ) 


CALL KERCAL 

IF (NS5 .EG. 0 ) GO TO 104 
WRITE (NWOUT* 7) 

DO 106 I=KOfKI" 

WRI TE < NWGUT * 5 ) ( J f Q < I f J ) f J-'KO * KF ) 
DO 115 K=2fKKF 
DXS=DX(K)+DX<K-1) 


D 1 < K ) =DX < K ) / ( DX ( K- 1 ) *DXS ) 

D2 ( K ) - ( DX ( K ) - DX ( K- 1 ) ) / ( OX ( K ) #DX < Tv~ 1 ) ) 


D3 ( K ) -DX ( K-l ) / ( DX ( K ) *DXS ) 

CONTINUE 

DO 510 K-l * KF 

FS ( K ) ”0 , 0 

DO 610 K-l f KF 

AXK=ABS < X ( K ) ) 

S X K = S 0 R T ( A X K * 1 2 - 1 . 0 ) 
l-IZ ( K ) -0 , 0 

IF (AXK ,LE. 1*0) GO TO 610 

HI ( K ) -ADS ( AXK*SXK ) - ALOG ( AXK+SXK ) 

CONTINUE 


CALCULATIE INITIAL MOD 


HERR -2 . 374*UB**0 ♦ 125/SORT ( 48 * *POE ) /C , 8 

DO 1990 IE=1fITE 

MTJ=2 

IF (IE * EG * 1) MTI-1 
IF (IE ,EQ, 1 ) IT OB- HERR 
DO 665 MT-MTI f MTMAX 
FMT=MT 

W T = ( F M T - 1 . 0 > * D W T + 1*57 0 7 9 6 3 


START II ME LOOP 


SINWT-SIN ( WT ) 

SINWT-ABS ( SINWT ) 

ITfi-1 

H0B2-H0B 

H0B1=H0B*0,95 

KF1-1 

I T OB -HO D1 


720 

C 


730 

C 

CALCULATE P(l) TO P(KA) 

740 

C 


74:1. 

1 22 

KFM= : KF~1 

750 


DO 123 K-l y K'FM 

760 

123 

PSA VE ( K ) =F < K > 

761 


CALL HCAL<1) 

766 


CALL MCAL ( 2 ) 

771 


IF ( NS1 * EG ♦ 0)00 TO 110 

7 76 


WRITE(NW0UTy6)H0B 

781 


WR I TE ( NtJOUT y 5 ) < K y H ( K ) y K = 1 1 KF ) 

786 

110 

CONTINUE 

791 


DO 120 K”1 y KF 

796 

120 

HS ( K ) =H ( K ) -HD ( K ) 

800 


KKO-KO-l 

805 


DO 531 K=lfKA 

810 


HH= ( H < K ) +H ( Kf 1 ) ) #0 , 5 

815 


TG(K)=1.0/HH**3 

820 


I'R ( K ) •- ( HH- 1 . )*TG(K> 

830 

581 

CONTINUE 

835 


IF(MT .EG. 1) GO TO 585 

840 


HDH” MODI.../ HOB 

845 


"DO 582 K* — 1 y K 0 

850 

582 

Ti l ( K ) = ( H ( K ) - HDH3KHL ( K ) ) 

870 


FS < KO ) -() . 0 

875 


DO 583 K=lrKKO 

880 


KK--KKO-K + 1 

885 

583 

FS < KK ) -0 . 5* ( TH < KK ) +TH < KK+ 1 ) ) *DX ( KK ) +F3 ( KK+ 1 ) 

890 

585 

CONTINUE 

891 


P ( 1 ) ~0 , 0 

892 


DO 586 l...”2 y KA 

893 


P<L)=P(L.-1 ) + <SINWT*FR<L-l )-< 1 .0/D7 ) * ( FS ( L- 1 ) +FS (L 

894 


.1 *T0 ( L - 1 ) ) *DX ( L- 1 ) /HOD **2 

895 

586 

CONTINUE 

896 

124 

I T ™ 1 

897 


IFCIIM » 0 T ♦ 1) 00 TO 590 

898 


CKI=0*<) 

899 


CK0=0 ♦ 0 

900 

590 

CONTINUE 

905 


IF (MI ,E0. 1) GO TO 372 

906 


HDH-HOBL/HOB 

907 


DO 371 K^KOy KI- 

908 

371 

TH ( K ) = ( I I ( K ) -HDH. t HI.. ( K ) ) 

909 

372 

CONTINUE 

910 


NNN--KF-KO 

911 


DO 340 K-lyNNN 

912 


KK=KF-K 

913 


TEMPI := 0.0 

914 


IF ( MT * EQ « 1) GO TO 37 3 

920 


TEMP' 1 --SUMS ( KK y TH y 1 ) 

922 

373 

T E M P = < H ( K K > - 1 , 0 > * S I N U T - T E M P 1 / D T 

924 


IF (TEMP » 0 T * 0.0) GO TO 340 

926 


W R I T E ( N U 0 U T 1 1 8 ) K K . H ( K K ) y T E M P 1 r T E M P 

928 


KFF=KK+1 

930 


GO TO 345 

932 

340 

CONTINUE 

940 

345 

CONTINUE 

941 


DO 417 K=KFF»KF 

942 

4 1 7 

P ( K ) =0 . 0 


1030 

C 


1040 

C 

CALCULATE FR ( K ) AND FS(K> 

1050 

C 


1055 

125 

CONTINUE 

1060 

126 

KKO-KQ-l 

1070 


DO 130 K“ 1 y KO 

1075 


HH=(H(K)+H<K+1 >>*0,5 

1080 


TG ( K ) -• 1 » 0/HH**3 

1090 


FR ( K > ~ ( HH--1 , ) *TG ( K > 

1097 

130 

CONTINUE 

1100 


IF(MT , EQ, 1) GO TO 140 

1110 


HDH-HOBL/HOB 

1120 


DO 135 K=1 .• KFF 

1121 

135 

TH < K > = ( I I ( K > -HDH*HL ( K > > 

1131 


IF ( l'TM » LE , ITMP > GOTO 560 

1132 


IF(MT , NE , MTP > GOTO 560 

1133 


WRITE < NWOUT r 16) 

1137 


WR I TE ( NWOUT t 5 > ( K r TH ( K > > K= 1 r KFF > 

1140 

560 

FS ( KO > -0 , 0 

1150 


DO 136 K-l f KKO 

1160 


KK--KKO-KM 

1170 

136 

FS ( KK >"0,5*(TH( KK > + TH < KK+ 1 > ) *DX ( KK > +FS ( KK+ 1 > 

1180 

140 

CONTINUE 

1190 


IF ( NS4 , EQ , 0 ) GO TO 139 

1191 


IF ( MT , EQ , 1) GO TO 139 

1200 


WR I TE ( NWOUT r 5 > ( K , FR ( K > t K = 1 r KO > 

1210 


WR I TE ( NWOUT r 5 > ( K t FS ( K > t K~ 1 , KO > 

1220 

139 

CONTINUE 

1340 


IF ( NS1 , EQ , 0 > GO TO 143 

1350 


WRITE ( NWOUT r 9 > IT r CKI r CKO 

1351 


IF ( NS6 ♦ EQ, 0) GO TO 143 

1360 


WR I TE ( NWOUT y 5 > <KfPS(K> r K~1 r KFF > 

1370 


WR I TE < NWOUT f 5) (K'rH(K) » K=1 r KFF ) 

1380 

143 

CONTINUE 

1390 


NN-KFF-KA 

1400 

C 


1410 

C 

A (.1 , 1 > y A ( 1 r 2 > 

1420 

C 


1490 

C 


1500 

C 

C( 1 ) 

1510 

C 


1520 

C 


1530 

C 

C(N> »N=1 >NN 

1540 

C 


1545 


KFM =KFF-1 

1560 


DO 170 K~1 r KFM 

1570 


HP ( K > = ( H < K > +H < K+l > > *0 . 5 

1580 


HP2 ( K > =HP ( K > **2 

1590 


HP3(K)=HP2(K>*HP(K> 

1595 


DP < K ) = < P ( Kf 1 > -P (K) > /DX C K > 

1605 

170 

CONTINUE 

1610 


DO .175 N~1 y NN 

1620 


K“N+KA“1 

1621 


C ( N > =--H0B**2*HP3 ( K > *DP ( K > + SINWT* ( HP ( K > ~ 1 , 0 > 



1622 


I F ( K . LT. K0) C(N)=C(N) -CKI 

1623 


IF ( K ♦ GE ♦ KO) C(N)“C(N)~CKO 

1637 

561 

CONTINUE 

1640 


IF<MT .EQ. 1) GO TO 175 

1641 


SUMTK ( K ) “SUMS ( K r TH y 0 ) 

1650 


C(N)=C(N)-1 .0/DT*SUMTK(K) 

1670 

175 

CONTINUE 

1680 

C 


1690 

C 

A(N»M) »N=1 »NNf M=1 f NN 

1700 

C 


1710 


Tl=l .5*C7*H0B/PI 

1720 


T2=H0B**2 

1730 


T3=C7/PI/H0B 

1740 


DO 200 N~1 y NN 

1741 


K=N+KA-1 

1742 


A(N»1)=0.0 

1743 


A ( N r NN ) =0 . 0 

1744 


I F ( K . L T . KO) A ( N y 1 ) == 1 » 0 

1745 


I F ( K * GE ♦ KO) A ( N y NN ) = 1 * 0 

1759 


NNM=NN-1 

1760 


DO 200 M=2yNNM 

1770 


J=M+KA-i 

1772 


I F < J ♦ GE ♦ KO) J=J+1 

1790 


DO 505 KG : ~1 y KFF 

1800 

505 

QQ ( KQ ) =0 ( KO y J ) 

1860 


A ( N y M ) =~HP2 ( K ) *DF : ' ( K ) * < QQ < K+ 1 ) +00 ( K ) ) *T 1 

1870 


I F ( J . EG . K+ 1 ) A ( N y M ) “A ( N y M ) +T2*HP3 ( K ) /DX ( K ) 

1880 


IF ( J . EQ . K ) A ( N y M ) “A ( N y M ) ~T2*HP3 < K ) /DX ( K ) 

1890 


A ( N y M ) =A ( N y M ) +T3*SINWT*0 ♦ 5* ( 00 < K ) +00 ( K+l ) ) 

1891 


IF ( NS7 .EG, 0) GO TO 562 

1895 


WRITE (NUOUTy 17) N y M y A ( N y M ) 

1900 

562 

IF (MT * EO . 1 ) GO TO 200 

1910 


A ( N y M ) = A ( N r M ) -T3*SUMS < K y 00 » 0 ) /DT 

1911 


IF ( NS7 ♦ EO * 0) GO TO 200 

1920 


UR I TE ( NWOUT y 1 7 ) N y M y A ( N y M ) 

1930 

200 

CONTINUE 

1940 

201 

IF ( NS2 ♦ EQ * 0 ) GO TO 210 

1950 


DO 205 KK==1 y NN 

1960 

205 

UR I TE ( NWOUT y 5 ) ( J J y A ( KK y J J ) r J J=1 y NN ) 

1970 


WRITE (NWOUT » 14) 

1980 


WRITE ( NWOUT y 5 ) ( KK y C ( KK ) r KK=1 r NN ) 

1990 

210 

CONTINUE 

2000 


l'DGT-'O 

2010 


CALL LEQT1F ( A y 1 y NN y 60 y C y IDGT y WORK r IER ) 

2020 


IF(NS1 .EQ. 0 ) GO TO 220 

2021 


NNM-NN-1 

2030 


WRI TE( NWOUT y 12) 

2031 


WRITE(NWOUTyll) C(l)y C(NN) 

2040 


WRITE ( NWOUT y 5 ) ( KK y C ( KK ) y KK=2 y NNM ) 

2050 

220 

CONTINUE 

2060 


CV==1 .0 

2070 


CKI”CK.I+C ( 1 ) 

2071 


CKO--CKO+C ( NN ) 

2072 


NNM=NN-1 


2090 


DG225 N-2 ? NNM 

2091 


K-NiKA™ .1 

2092 


IF(K .GE. K0) K=K+1 

2100 

2120 


DPS ( K ) =C ( N ) 

2121 


IF ( DPS ( K ) ♦ LT * -3» 0)80 TO 499 

2130 


I F ( ABS < DPS ( K ) ) -EPSP ) 222 , 222 ,221 

2140 

22 L 

GV=0.0 

2150 

9 99 

PS ( K ) "PS ( K ) +DPS ( K ) 

2160 


P < K > "PD ( K > f PS < K ) 

2180 

n n k:; 

CONTINUE 

9 9 9 C\ 

■V urn Am A** \X 

c 


2230 

c 

CHECK FOR N-R CONVERGENCE 

2240 

c 

2250 


■IF (CO ♦ E Q , 1.0) GO TO 240 

2260 


IF (IT .GT. I TP) GO TO 500 

2270 

c 


2280 

c 

N-R NEXT ITERATION 

2290 

c 


2300 


IT-IT'M 

2310 


CALL HCAL ( 1 > 

2311 


NS 2=0 

2312 


NS4=0 

2315 


GO TO 125 

2320 

240 

CONTINUE 

2340 

C 


2350 

C 

N-R CONVERGED , NEXT MAIN INLET ITERATION 

2360 

C 


2390 


PH I -OKI -CKO 

2400 


GO TO ( 594 , 595 ) , KF1 

2410 

594 

F'F.1. -PHI 

2411 


W R I T E ( N W 0 U T , 2 1 ) I TM ? HOB, P H I 

2420 


HOB- HOB 2 

2430 


KF1 -2 

2450 


ITMTTMil 

2460 


GO ro 122 

2470 

595 

PF2-PHI 

2490 


D P H I = ( p F 2 - P F 1 ) /(HO B 2 - H 0 B 1 ) 

2500 


DHOB=- PHI/DPHI 

2510 


WRITE ( NWOUT, 10) ITM, HOB, PHI, DPHI, DHOB* IT 

2520 


HOB-HOBiDHOB 

2530 


IF (ABS (PI- 1 1 ) . LT, EPSH) GO TO 460 

2540 


I F ( I TM .GT. ITH) GO TO 495 

2541 


IT M- IT Mi 1 

2542 


H0B1 -HOB 2 

2543 


PF1=PF2 

2544 


H0B2-H0B 

2545 


GO TO 122 

3350 

460 

CONTINUE 

3351 


W R I TE(NU, 8 )MT,WT, H 0 B 

3950 


WR I TE ( NW0UT , 8 ) MT , WT , HOB 

3951 


NWG-NWOUT 

3952 


IF (MT .EQ. MTMID) NWQ=NWOUT 

3953 


IF(MT ,E0, MTMAX) NWO=NWOLJT 


3960 

3961 

3962 
3965 
3266 
3967 
3° 68 

3971 

3972 

■7 Q 7 3 

3974 

3990 

4000 

4010 

4020 

4030 

4031 

4032 

4033 

4042 

4043 
A / j /X 

4031 

4052 

4053 

4054 

4055 

4056 

4060 

4061 

4062 
4 0 6 3 
4064 
4 0 6 3 
4 0 6 6 

4067 

4068 

4069 

4070 

4071 

4072 

4073 

4074 

4075 

4076 

4077 

4078 
4 079 
4080 

4104 

4105 
4 .! 06 
4 ! 07 
4108 
4 I OS ’ 
41 1.0 
4 1 1 . .1 
4 1. 1 2 


WR I TE ( NWR r 5 ) ( K t P < K ) r K=1 y KF ) 

WR I TE < NWG » 5 ) < K » H < K ) y K= 1 s- KF ) 

SUMKFF ( MT ) =SUMTK ( KFM ) 

HOBST ( MT ) =HOB 
DO 465 K== 1 y KF 
PS TOR < MT r K ) =P ( K ) 

465 HSTOR ( MT ? K ) =H ( K ) 

IF(NS8 .EQ. 0 ) GOTO 640 

WR I TE ( NWOUT y 5 ) ( K > TH ( K ) y K" 1 y KFF ) 

WR I TE ( NWOUT 1 5 ) ( K y SUMTK ( K ) y K-KA y KFM ) 

640 CONTINUE 

CALI... HCAL < 2 ) 

DO 270 K=1»KF 
l-IB ( K ) “H ( K ) HD < K ) 

1 1 SI... ( K ) =HS ( K ) 

HI... < K ) =H ( K ) 

270 CONTINUE 
HOBL. •"•TIOB 
GO TO 121 

495 IF ( KF5 .EQ. 1) GO TO 500 
WRITE (NWy 26) 

NWOUT "6 

NS1=1 *NS2=0 $ NS3=1 *NS4= 1 *NS6=1 
KF5 = 1 
GO TO 122 

499 3 F' ( K F 5 ,EQ. 1) GO TO 500 

WRITE (NWy 27) 

NWOUT-- 6 
NS I = I. 

N B 4 1 
KF.3=1 
00 TO 125 

121 IF ( MT . GT . 1) GO TO 655 
H IN-HOB 
DO 650 K--lyKF 
650 1 1 N T ( K ) :::: I I ( !\ ) 

655 I F ( M T , NE , MTMAX ) GO TO 660 
HI.. AST- HOB 
PI-1 1 A=H IN-HI... AST 

660 IF ( MT , NE ♦ MTMID ) GO TO 665 
DO 66* K= 1 y KF 
PTEM(K)=P(KF+1-K) 

HTEM(K)=H(KF+1-K) 

666 CONTINUE 

DO 667 K : -l y KF 
P ( K ) =FTEM ( K ) 

667 I II.. ( K ) ;::: l I TEM ( K ) 

WR I TE ( NWOUT y 5 ) ( K y HL.. ( K ) y K- 1 y KF ) 

665 CONTINUE 

DO 680 K~ : 1 ? KF 
680 HL(K)=H(K) 

IF (IE , E Q . 1) GO TO 670 
QF2~PHI A 

D P H I A ( Q F 2 - Q F 1 ) / ( H I N 2 - MINI) 

DHIN— PHIA/DPHIA 
WRITE (NWy 22) 

WR I TE < NW y 23 ) I E y H I N 1 r H I N2 y PH I A y DPFI I A r DH I N 
QF1--QF2 


4113 

4114 

4115 

4116 

4117 


MINI HIN 2 
HIN=HIN-FDHIN 
HIN2=HIN 
IF(ABS<DFHN> 


HOBL-ITIN 


* L.. T ♦ 


EPSE) GO TO 2000 


4118 

4119 

4120 

4121 

4122 

4123 

4124 

4125 

4126 

4127 

4163 

4164 

4165 

4166 

4167 

4168 

4169 

4170 

4171 

4172 

4175 

4176 

4177 

4178 


HOB-HIN 
GO TO 199() 

67() HIN1=HIN 

H I N 2 IT IN/:. * 2 5 
QF1=PHIA 

WR I TE ( NW y 24 ) I E y H I N y PIT I A 
HIN-HIN2 
HOBL=HIN 
HOB=HIN 
1990 CONTINUE 
500 CONTINUE 
2000 CONTINUE 

SUMKFF ( 1 ) ~0 . 0 

BO 2010 MT-.1. yMTMAX 

FMT--MT 

WT“(FMT--1 . >*DWTM .5707963 

WRITE < NW y 25 ) MT'y WT > HOBST ( MT ) y SUMKFF (MT) 
WRITE(NWyS) <Ky PSTOR < MT y K ) y K™ 1 y KF ) 

WR ITE(NUyfv) ( K y HSTOR < MT » K > t K= J y KF ) 

2010 CONTINUE 
STOP 

1 FORMAT (72H 

2 FORMAT <1615) * 


FORMAT' ( 7E10 . 3 ) 

F ORMA r </4X 1 4HP0E= t El 2 . 5 y 5X y 3HUB-- y El 2 . 5 y 5X y 4HDWT= y El 2 ♦ 5/ ) 

F ORMAT ( 8 ( 1 X y I 2 ? E 12.5) ) 

F ORMA T ( /4X r 4 ITU ( K ) y 5X. 5 H IT D ( K ) y 4 IT U O B ;■= y F12.5/1 

F0RMAT(/4X»6HQ(Ir J)/) 

I- ORMA I ( /4X y 3MMT= y 15 y 5X y 3HWT“ y El 2 . 5 y 5X y 4HH0B“ y E.1.2 . 5/ ) 

FORMA I ( / 4Xy 3HIT= y 15 y 4Xy 4HCK I “= y El 2 . 5 y 4Xy 4 IT CKO"' y E12.5/) 
F ORMA I < /4X y 4FIITM- y 15 y 5Xy 4HH0B= t E12 . 5 y 5X y 4ITPHI- y E.1.2 ♦ 5 y 5X y 


1 5HHPITP 


5FIDIT0B= y E 1 2 ♦ 5 y 3FI 1 1 


r~ w * w. , ’ * — ■ • • « w r i... * * w r wl IJ, I “ X.C. / 

FORMAT <S(lXyE12.5) ) 

2 FORMAT ( 4X y 16HDCKI * BCKO , DPS (K) ^ 

FORMAT < 2X y 5HC < KK ) ) 

FORMAT </5X»4HKFF=y I5/> 

FORMAT ( 2X y 5HTH < K ) > 

F ORMA T < 7 < 1 X r 2 1 2 y E 12.5)1 
F ORMAT ( 15 r 3E.1.2 ♦ 5 ) 

1HNy AX * IHKr 4X y 1 IT J y 4Xy 6H0<KyJ)y 

1 7Xy 8H0 (K + l yj) y 5X y 6HA(NyM) ) 

FORMAT ( 415 y 6<lXyE12.5)) 

FORMAT < /4X y 4FI ITM“y I5y5Xy 4HH0B= y El 2 . 5 y 5X y 4HPHI = y El 2 . 5 ) 

~ ' 2HIE ' 3X ' 4HHIN1 y 9X y 4HHIN2 y 9X y 4HPH1A y 9X y 5HDPHIA y 

1 oa y 4HLIHIN ) 

5 FORMAT <4XyI5y5(lXyE12.5)/) 

\ FORMAT ( /4X y 3HIE= y 15 y 3FIHIN y E12 . 5 y 4HPHI A y El 2 . 5/ i 

3HWT =’ E12.5,SX, 4HH0IW, 

1 El a-.. ♦ %j f /HSUMKF* F- ?E12»5/) 

> FORMAT ( 13HITMEXCEEDEB /) 

’ FORMAT ( 13HDPS NEGATIVE /) 

END 


4 2 9 2 S U B R 0 U T I N E K E R CAL 

4293 COMMON HOB y HOB A (25) r WT r DT r KA r KO r KF t KR r KFF r UB r POE 

4310 COMMON X ( 80 ) y DX ( 80 ) r Q < 80 y 80 ) y 0 JA ( 80 ) 

4320 COMMON H ( 80 ) « HD ( 80 ) r HS < 80 ) r HL (80) y HSL ( 80 ) 

4330 COMMON P ( 80 ) y PD ( 80 ) , PS ( 80 ) t DP ( 80 ) t DPS ( 80 ) » PSAVE ( 80 ) 

4340 COMMON FR ( 80 ) y FS ( 80 ) y D1 ( 80 ) y D2 ( 80 ) y D3 ( 80 ) 

4350 COMMON TA ( 80 ) y TB ( 80 ) y TG ( 80 ) y ID ( 80 ) y TE ( 80 ) y IF ( 80 ) y TO ( 80 ) 

4360 DO 1 1=1 y KF 

4370 DO 1 J=l, KF 

4380 1 G(IyJ)=0.0 

4390 KKF=KF-2 

4400 BO 8 K=1 y KF 

4410 G ( K y 1 ) 0 ♦ 0 

4420 F5=X(K) 

4430 DO 8 J = 1 y KKFy 2 

4440 U=X (J)-F5 

4450 U2--X ( J f 2 ) -■ F 5 

4460 AU=ABS(U) 

4470 AU2= ABS(U2) 

4480 IF(AU) 50 y 51 y 50 

4490 50 AU=ALOG(AU) 

4500 51 IF ( AU2 ) 52y 6y 52 

4510 52 A U 2 = A L 0 0 ( A U 2 ) 

4520 6 DJ=X( J+l)-X( J) 

4530 F2=3.0*DJ 

4540 UQ=U*U 

4550 U2Q==U2*U2 

4560 FK=UQ* ( Al.J-1 ♦ 5 ) *0 . 5 

4570 FK2=U2Q* ( AU2- 1 * 5 ) *0 * 5 

4580 FKB=U* ( FK-UQ/6 ♦ 0 ) -U2* ( FK2-U2Q/6 * 0 ) 

4590 Q ( K y J ) = ( ( -3 . 0*FK--FK2 ) /2 ♦ 0-FKB/F2 ) /D J-U# ( AU- 1 . 0 ) + Q ( K y J ) 

4600 G ( K y J+l ) = ( 2 ♦ OK ( FK+FK2 ) +2 ♦ 0KFKB/F2 ) /D J 

46 1 0 G ( K y J+2 ) = ( ( - FK-3 , 0KFK2 )/2. 0-FKB/F2 ) / D J+U2* ( AIJ2- 1 . 0 ) 

4620 8 CONTINUE 

4630 KFP=KF+1 

4640 DO 100 J=lyKF 

4 6 5 0 1 0 0 G ( K F P y J ) = G ( K R y J ) 

4660 DO 300 K=1 y KF 

4670 DO 300 J=l» KF 

4680 300 G(Ky J)=Q(K» J)~G(KFPy J) 

4690 RETURN 

4700 END 


4710 SUBROUTINE HCAL (KFORK) 

4720 COMMON HOB y HOB A ( 25 > r WT r DT r KA r KO y KF y KR t KFF y UB r POE 

4730 COMMON X < 80 ) , DX ( 80 ) y 0 ( 80 , 80 ) y Q JA ( 80 ) 

4740 COMMON H ( 80 ) - HD (80) , HS ( 80 ) , HL ( 80 ) y HSL ( 80 ) 

4750 COMMON P ( 80 ) y PD ( 80 ) , PS (80) y DP < 80 ) y DPS < 80 ) y FSAVE ( SO ) 

4760 COMMON FR ( 80 ) ?FS(80) t D 1 ( SO ) , D2 (BO) r D3 ( 80 ) 

4770 COMMON TA ( SO ) y TB ( 80 ) y TC ( 80 ) y TD ( 80 ) y TE < 80 ) y I F ( 80 ) y TO ( 80 ) y TH ( 80 ) 

4771 COMMON GQ ( 80 ) y HZ <80 > 

4780 PI -3. 141593 

4790 07“ 1 6 * 0*PQE**2/SQRT ( 48 , *UB ) 

4800 C8“C 7/HOB 

4810 IF (KFORK , EQ ♦ 2) GOTO 20 

4820 DO 10 K=1 y KF 

4821 IF ( ABS ( X ( K ) ) .LE, 1 ♦00001) GO TO 5 

4830 H<K)=0. 

4840 DO 1 J=lyKF 

4850 1 H ( K ) ■•■H ( K ) -H P ( J ) -PD ( J ) ) *Q ( K y J ) 

4860 H ( K ) = 1 . +C8* ( 0 ♦ 5*HZ < K ) -H ( K ) /P I ) 

4361 GO TO 10 

4862 5 H<K)=0.0 

4863 DO 6 J=lyKF 

4864 6 I I ( K ) = H ( K ) + ( P ( ...I ) - P D ( J > ) * Q ( K y J ) 

4865 H ( K ) = 1. . 0-C8*H ( K ) /P I 

4866 10 CONTINUE 

4870 GO TO 30 

4880 20 DO 25 K~1 y KF 

4390 HD ( K ) - : 0 ♦ 

4900 DO 21 J=iyKF 

4910 21 HD ( K ) “HD < K ) +F : 'D ( J ) *0 ( K y J ) 

4920 25 HD ( K ) =1 ♦ O+CB# ( 0 ♦ 5*X ( K ) **2“HD ( K ) /PI ) 

4930 30 RETURN 

4940 END 


4950 

4960 

4970 

4980 


FUNCTI ON SLJMN < K 1 y K2 y A ) 

DIMENSION A (80) 

COMMON HOB y HOBA ( 25) y WT y DT y KA y KO y KF y KR y KFF y UB y POE 
COMMON X ( 80 ) y DX ( 80 ) y Q ( 80 y 80 ) y OJA ( 80 ) 


4990 KK2-K2-1 


5000 

50.1.0 

5020 10 

5030 

5040 

5050 

5060 


SUM-0 ♦ 0 

DO 10 K-Kl y KK2 

SUM-SIJMTDX ( K ) #A ( K ) 

SUMN=SUM 

RETURN 

END 


5070 


FUNCTION SUMS ( K1 y A y NEXIT ) 

5080 


DIMENSION A (80) 

5090 


COMMON HOB r HOBA ( 25 ) y WT y DT y KA y KO y KF y KR y KFF 

5100 


COMMON X ( 80 ) y DX ( 80 ) y 0 ( 80 y 80 ) y Q JA ( 80 ) 

5110 


KK0=K0-1 

5120 


KK1-K.1. f 1 

5130 


KK2=K1-1 

5140 


I F ( K 1 ,LT. KKO ) GO TO 10 

5150 


I F ( K .1. * EQ . KKO) GO TO 20 

5160 


I F ( K i .EG. KO) GO TO 30 

5170 


I F ( K 1 .GT. KO) GO TO 40 

5.1.80 


SUMS-0 . 0 

5.1.90 


RETURN 

5200 

.1.0 

SUM=DX ( K 1 ) /4 . 0* ( 0 . 5# A ( K 1 ) + 1 . 5# A ( K 1 + 1 ) ) 

5210 


DO 100 K-KK1 y KKO 

5220 

100 

SUM-- SUM1DX ( K ) /2 . 0* ( A ( K ) +A ( K+ 1 ) ) 

5230 


SUMS=SUM 

5240 


RETURN 

5250 

20 

SUMS=DX ( K 1 ) /4 . 0* ( 0 . 5* A ( K 1 ) + 1 . 5* A ( K 1 + 1 ) ) 

5260 


RETURN 

5270 

30 

SUMS--DX ( KO ) /4 , 0* ( 1 . 5*A ( KO ) +0 . 5*A ( K0+ 1 ) ) 

5280 


RETURN 

5290 

40 

SUM-0 . 0 

5300 


DO 200 K-KO y KK2 

5310 

200 

SUM-SUM-DX ( K ) /2 . 0 * ( A ( K > +A ( K+ 1 ) ) 

5315 


IF (NEXIT .EG. 1) GO TO 210 

5320 


SUM-SUM-DX ( K ) /4 . 0* ( 1 . 5* A ( K1 ) +0 . 5*A ( Kl+1 ) ) 

5321 

210 

SUMS -SUM 

5330 


RETURN 

5340 


END 


Input Instructions 


Card 1 Format (72H • • • ) 

Identification card for the run 

Card 2 Format (1615) 

NRUN — Number of runs to be made 

Card 3 Format(16I5) 

KA — k^, the grid which divides the contact into two regions. 

For k < k^, the pressure is determined by integrating the 
Reynold's equation directly [ 13] . For k > k , the Newton- 
Raphson method is used. 

K<j> — k Q , the grid number at the contact center. 

KF — k^., the last grid. 

KR — k , the grid at which H = 1. 

K 

KFF — k ff , the grid at which the film terminates. 

NHERR - Not used in RODSLJ, set NHERR = 1. 

Card 4 Format(16I5) 

NS1 to NS8 — diagnostic controls, set NS1 to NS8 = 0. 

Card 5 Format (1615) 

ITH — Maximum iterations allowed for ITM. Recommended ITH = 15. 

I TP — Maximum iterations allowed for IT. Recommended ITP = 25. 

ITE — Maximum iterations allowed for IE. Recommended ITE = 4. 

MTMAX — (MTMAX-1) represents the total number of time intervals 
within half of the reciprocating cycle. 

MTMID — The time step MT corresponding to mt=ir, or the end- 
stroke position. 

Card 6 Format (7E10 .3) 

EPSH — Maximum error allowed for 4>__ ,(<{>__= I c *' — c'l). Recommended 

H H I O 

EPSH = 0.0005. 

EPSP — Maximum error allowed for AP, . Recommended EPSP = 0.0005. 

k 

EPS E — Maximum error allowed for <J> E , (<|> E =|h o ^ - h ^ ^ | ) . 
Recommended EPSE = 0.005. 


X( 1) — The x coordinate for the first grid. 


Card 7 Format (7E10. 3) ___ 

DX(K) (K=l , (KF-1) ) A set of cards for Ax k = x k+1 - x k * 


Card 8 Format (7E10. 3) 
P(f>E - p /E' 


DWT — Amt (radians) 

FL4B - S/4b 

Card 9 Format (7E10. 3) 

P(K)(K=1,KF)- A set of cards for the intitial guess of pressure 
profile P k » 


Sample Input Data 


The following is a set of sample input data for run 51 listed in 


Table 2. 


10 ROD A 

30 




20 1 





30 3 

21 

4.1 2.1 

32 

.1 

40 0 

0 

0 0 

0 

0 

50 15 

25 

4 19 

10 


60 

.0005 

.0005 


.005 

70 

0.2 

0.2 


0 . 03 

80 

0.01 

0.01 


0.01 

90 

0.06 

0.06 


0.2 

100 

0.2 

0.2 


0.2 

110 

0.01 

0.01 


0.01 

120 

0.01 

0.03 


0.03 

130 

0.10 

5 . 00 E -07 

. 1745327 

140 

0.0 

0.0 


0.0 

150 

.035 

.09 


.1410 

160 

♦ 39.19 

.5102 


0.6 

170 

.9797 

.9165 


0.8 

180 

0 . 275 

0 . 235 


0.18 

190 

o 

* 

o 

0.0 


0.0 


0 0 


. ♦ 5 




0.03 

0.01 

0.01 

0.01 

0.01 

0.01 

0 . 02 

0 . 02 

0.2 

0.2 

0.2 

0 . 2 

0.06 

0.06 

0.02 

0 . 02 

0.01 

o 

* 

o 

0.01 

0.01 

0.2 

0.2 



5.00 




0.0 

0.01 

.02 

.035 

.200 

.2431 

.28 

.3411 

0.8 

.9165 

.9797 

1.0 

0.6 

.5102 

.3919 

.3411 

0.10 

0 . 035 

0.0 

0.0 

o 

o 

o 

o 

0.0 



Sample Output Data 


A sample of output data based on the input data for Run 51 is included. 
The first part lists the values of h Q for different stroke positions 
for each iteration IE. For this sample run, the IE loop converges 
after three iterations. Thus, results of h Q for IE = 3 are taken to 

be the final solution for Run 51. In the second part of the output, 
there are two arrays of data listed under each stroke position uit. 

The first array is the pressure profile, and it is followed by the 
film thickness profile. 

The value listed as SUMKFF represents SUMS (KFF , TH , 0) which can be used 
in the future for calculating the leakage rate at each stroke position. 



UT= 

.15708E+01 

HOT" 

, \ •. r:;r too 

m r 

WT= 

• 17453E+01 

MM): 

• loo 

'!• 

WT= 

• 19199E+01 

t ' f n : 

•:> -or i-oo 


WT= 

.20944E+01 

HUP • 

4 :■! ;'5AOPT-fO'') 


WT= 

.22689E+01 

HO (3 

r oo 

.‘if 

UT= 

• 24435E+01 

HOB. 

, 19579F+00 

Ml 

UT= 

• 26180E+01 

HOB- 

- 1 O'! 6? Ft- 00 

r-rv 

WT= 

.27925E+01 

HOB- 

. 17439E+00 

M-: 

R 

h- 

3 

I 

•29671E+01 

non- 

. 16704E+00 

MT ■ * ■- 

WT= 

. 3141&E+01 

HOB-- 

. 16331E+00 

iS j . • 1 ! 

WT= 

.33161E+01 

HOB- 

. 1.61 18E+00 

HI 

WT= 

• 34907E+01 

HOB'-- 

, J 5442E1 00 

• ; 1 • : 

WT= 

.36652E+01 

’ HOB ■ 

. 14J.73E+00 

MI ■ i ' 

WT= 

•38397E+01 

HOB • 

- 139A7E+00 

1 1 1 : 

WT= 

• 40143E+01 

HOB-- 

. 1. 0-1 4 6c -f 00 

• i _ 

UT= 

.418B8E+01 

HOB- 

1 760 IE 4-00 

■ 

WT= 

• 43633E+01 

HOB 

, | v,iPf ,|-.}00 

til 

WT= 

• 45379E+01 

HOB ■ 

» BOS 1 K 1 00 

! - 

«!=■ 

• 47124E+01 

Hub 

• ’ 1 1 c !• 0 0 


1 - ■ 


( 


■I 



1 


• I'l 


a 


UT= 

• 17453E+01 

MOB- 

, :: l 9O2E + 00 

UT= 

. 19199E101 

HOP 

. 1. f)?8E 100 

UT= 

• 20944E101 

1 !■'>.' 

,21 3 4 BE. 100 

ur= 

.226B9E101 

HUB 

, 20H72E100 

UT= 

• 24435E101 

1 1 n • • 

. 1 05721: 100 

WT= 

• 26180E+01 

H('I.v ■ 

■ 1 0465E100 

WT= 

.27925E+01 

1 101' 

. 17436E100 

WT= 

•29671E+01 

I I! O' ^ 

, 1.6 7 07 El 00 

UT= 

•31416E+01 

hob 

, 16329E100 

UT= 

.33161E101 

HOB" 

, 161 1 A El 00 

UT= 

.34907E+01 

HOB - 

. 15440E100 

UT= 

• 36652E101 

HOB 

, 141 71 El 00 

WT= 

•38397E+01 

1 IOB" 

. 1 396 6E 100 

WT= 

• 40143E101 

1 !0B* 

. 1 'll -1 lit:. 100 

WT= 

.418B8E101 

HOB 

.1760 IE 100 

WT= 

• 43633E+01 

HOB 

. 19485E100 

WT= 

. 45379E101 

HOB 

, 20317E100 

UT™ 

> 47124E101 

HOB 

, 1.3? IE 100 


H : M2 ! 

, 2159:1. E + OO 

'!H. 

DP 11:10 BHTM 

i "0E--0' 7 , 99999E100 


UT= 

* 17453E+01 

HUP- 

♦ 19258E+00 

WT= 

• 19199E+01 

hop - 

. 20439E » 00 

UT = 

• 20944E+01 

HOB 

, :o .-'('.of t o(i 

UT= 

•22689E+01 

HOP - 

. : : :0'>:.6L too 

WT= 

> 24435E+01 

HOP 

. L9477E+00 

WT= 

• 26180E+01 

HOP 

- 1.84.1.2E+00 

WT= 

< 27925E+01 

|lt (P: 

. :l. 740:1.171-00 

UT= 

•29671E+01 

HOP 

. 16675E+00 

WT= 

•31416E+01 

HOB- 

. 16303E+00 

WT= 

• 33161E+01 

HOB- 

. 16087E+00 

WT= 

• 34907E+01 

HOP = 

. 15404E+ 00 

UT= 

.36652E+01 

HOP = 

. 14138E+00 

UT= 

•38397E+01 

HOP:.-: 

. 1 3?47P 400 

UT= 

•40143E+01 

HOP 

. I.5438E 400 

WT= 

•41888E+01 

HOi;- ■ 

I 7599E+00 

WT= 

» 43633E+01 

HO P: : 

- 1 9485E-400 

WT= 

• 45379E+01 

HOP 

-2001 7E-400 

WT“ 

• 47124E+01 

HOP 

< 21591E-400 


Lift).: 

pi i r. i 

0I-H.1 A 


3 •' . r-^/BOliT-fOO • . 1. »9??99E+00 ♦ 38.108!-. -01 


t 



1 


MT“ 

j. UT“ 

. 157087+01 

nop'- 

'J . 

ie+oosumkff- 

0. 









! •' . * . 

,, 

■ 5 . •? ( T r. 

3 

,17707i. (M 

■1 

• 3437317-01 

5 

.6386211-01 

6 

.7047611-01 

“7 

/ 

, 958727-01 

') 

. 1 1 62 2E 

100 


• In 

(00 

1 :l. 

, 19346F + 00 

12 

. 22295E+00 

13 

, 252897+00 

14 

.311507100 

15 

.364907+00 

1 6 

, i • '.\sr 

[■ 0 0 


■ I ; > 

’ 9 7 8 8 E: ( 00 

19 

.•vcirir. (oo 

20 

. 979627+00 

21 

. 100007+01 

2 9 

. 9002411( 00 

23 

. 9.1 70011 100 

24 

, 80 a '-'Of 

! *00 


, -A -’. 

; i i ?■'!-■( 00 

27 

. 3956 Of 1 00 

23 

.3 18 -16 El 00 

29 

. 2867311 + 00 

30 

. 2 433.21:. ( 00 

31 

. 18.1.547 + 00 

3 2 

. 10 5f: 

-0 1. 

:vs . .i ■■■ i 

■t.- ; v i 


35 

0. 

36 

0. 

37 

0. 

38 


39 

0, 

40 ( 

, 


A 1. ' ' - 

! > ‘ • 


’32487+02 

3 

,633!ioi:r (oi 

■1 

. 4 593 5 El 01 

5 

.322937+01 

6 

.20640E+01 

7 

. 254487+01 

8 

, 22’ 1411 

! 0 1 


• U ’ !• n * 1 • ■ 

01327:1 or 

11 

. 1 /»7V 11 (01 

12 

. 15524 El 01 

13 

.1451711101 

1.4 

. 131.1 IE +01 

15 

. 12260EI0 l 

16 

, :!. 1. 2 127 

+ «) (. 

•j. ■’ .1 

! ' ! ' 

o ; r } o i 

19 

. 102 19 51- 01 

20 

. 10106E+01 

21 

.100001.101 

'■> '■) 

, 9923 4 Ell- 00 

23 

, 98058 E 4 00 

.2 4 

.9700'i.F 

! o 0 



■.:> ! ••■TK'O 

77 

.9189617 + 00 

20 

. 70 093 II ( 00 

29 

. 879297+00 

30 

.8500811+00 

31 

. B0883E+ 00 

3 r.' 


( 00 


■ I ; v- 

' 55031 ::: (oo 

35 

, 1.2 '2 5?;. (0 1 

36 

. 156327+01 

37 

. 19603E+01 

38 

. 344547+01 

39 

.530177+01 

40 

. 2 46 ion 

+ 02 

MT - 

2 UT" 

. 17453E+01 

IIOB . 2 

90 

;-e l ooniJMKrn.: 


233057-01 











'.‘1*7077-02 

3 

. 1 7 4 9 8 E 01 

. 4 

. 33957E-01 

wJ 

. 633157-01 

6 

.7790617-01 

7 

. 95298E-01 

8 

. 1.15677 (00 

1 V-: ' 


, 16 4927(00 

1.1 

. 1 93 1 3 II +00 

1.2 

. 2227217+00 

13 

.252777+00 

14 

.31 1527100 

15 

.365087+00 

1.6 

. 497497 +00 


i im 

. 79793E + 00 

19 

.9151411+00 

20 

. 979647+00 

21 

. 100007+01 

6':> 

. 980197+00 

23 

.916937+00 

24 

, 800 80S +00 



.51 1. 6317 + 00 

27 

.3948317 + 00 

28 

. 3 4 761.7+00 

29 

.286917+00 

30 

. 244517+00 

31 

.183367+00 

32 

,*>89 247 -01. 

| ..•••• : 

t o i / .i 

, 17- 49 2 E- 01 

35 

0. 

36 

0. 

37 

0. 

38 

0. 

39 

0. 

40 v 

> .. 


, n 

- 256097+02 

3 

.6417517+01 

4 

. 464727+01 

5 

.326327+01 

6 

, 28929E+01 

7 

.256937+01 

8 

,729947+01 


' <r\ 10 

■ .1861.47+01 

11 

. 16939E + 01 

12 

. 156597+01 

13 

, 146457+01 

14 

. .13231. E+01 

15 

• 12376E+01 

1.6 

. 1.13557 (01. 

1 - v: ■ 

. : M. *01. i« 

, 104927+01 

19 

. 102837+01 

20 

. 10155E+01 

21 

.100007+01 

22 

•98824E+00 

23 

• 972267+00 

24 

. 956587. 1 00 


: o '>0 

.915497+00 

27 

.8989617 + 00 

28 

.88998E+00 

29 

.863327+00 

30 

. 835527+00 

31 

. 794 12E (00 

32 

,766577+00 


. r *. o •••; .a 

. 93546E+00 

35 

. 1209117+01 

36 

. 155557+01 

37 

.19592E+01 

38 

. 346787+01 

39 

. 535267+0.1. 

40 

. 2 4952F+02 


WT= . 19199E+01 HOB- . 2 1 8 2 87+008 I J M K F F = - 


.1 0, 


r 23767E -02 

3 

. 1673817-01 

4 

.328027-01 

5 


- 00 

t 16JA7EV0Q 

1.1 

. 1919311+00 

12 

,221807+00 

1.3 

S 7 ■ 7' : 

■ ; T'.j i p 

. .’978911+00 

19 

.9154 1.17+00 

20 

. 97962E+00 

21. 

■ 

i ->o 

. 51 1. 7 27: + 00 

27 

. 395257+00 

28 

.347967+00 

29 

, 1 

■■•'.vl •: ) - 

, 142417-01 

35 

0. 

36 

0. 

37 0 

\ l 0 i 

.1 +'■■•'' ' ■ 1 

? 

■ 25665E+02 

3 

.639077 +01 

4 

. 46.1 lOE+O.l 

5 


? - ■ ! 0 

: 18 1.5 .1.7. +01. 

11 

. 1648817+01 

12 

. 15222E+01 

13 

i. • i ! 

: - vi ! v 

1 9323E+01. 

.1.9 

,10.18.1! +0 1 

20 

. 10094(7101 

21 

,J . 


.92937E +00 

27 

. 9 13237 100 

28 

.903187+00 

29 

VI . p. 1 - 

■ :■ -o ,-v-i 

: 30 3E f 00 

35 

, 122467+01 

36 

. 157377+01 

37 


579537-02 


617147-01 

6 

.762027-01 

7 

.9354317-01 

8 , 

, 1 13947+00 

252127+00 

14 

, 31 1 3.1.EM00 

15 

.365097+00 

16 . 

, 1973 47 + 00 

1.000017+01. 

2 2 

.980217+00 

23 

.916957+00 

24 

,8008 .’7 + 00 

286387+00 

30 

.243137+00 

3.1 

.181357+00 

32 

, 976837 01 


38 

0, 

39 

0. 

40 0. 


321977+01 

6 

.284767+01 

7 

.252297+01 

8 . 

.224557+01 

142257+01 

14 

. .1.28467+01 

15 

. 120227+01 

16 . 

, 1. I 0657 + 01 

100007+01 

-j 

, 992367 + 00 

23 

.979887+00 

24 

,'?/>■■• 8 5F+ 00 

873717100 

30 

.844367+00 

31 

.802457+00 

32 

, 776057 + 00 

197867+01 

58 

.349237+01 

39 

.538357+01 

40 

.250 4 77 +07 


- 1 1 , 



f1T~ 4 WT= . 20944E+01 HOB" . 2l34f?E+00F.IIMKFF= 


t . 



,29037!~ 02 

3 


■ 15A56E-01 

4 

. 31065E-0 1 


, W. 

* i i 

, 1.61 39EM-00 

it 


. 1 9032E + 00 

12 

. 22070E+00 

1 

| . . I- •! • n . 

•; 2 

.79793E+00 

19 


. 9 .1 542E+00 

20 

• 979AAE+00 

' ' 

t 4 : ■ 


.31 1.9 OEM 00 

2 7 


. _? ?J 5A2E + 00 

28 

. 34828E+00 

•1 1 

: 9A45E -0 1 

0 1 - 

. J 240 IE- 01 

35 

0 

• 

36 

0 . 

' : A |-. j,'. ' 

•) 

. 26193EM02 

3 


. 6453 4 EM- 01 

4 

. 46293E+01 

V 

1 -'MME + Ot 

1 0 

. 17757E+01 

11 


. .1.6 092 EM 01 

12 

. 14836E+01. 

:! '■ 

; •! ; • i 

■ B 

, 10 184 EM 01. 

19 


. 10096EM01 

20 

. 10047EM01 

c: •'.> 

''•.''loot: +■"'■•) 

2 /> 

.941 7 3 EM 00 

27 


.9 2553 EM 00 

28 

.91445E+00 

1 ! 

■' 1 '.* *'> <r_> 

\ A 

. 9 A El 7 3 EM- 00 

35 


. 1 2475EM 0.1 

36 

. 1.6019E+01 


MT~ 5 

WT= 

• 22689E+01 


HOB- . 20 

'572E+00SUMKFF-- 

1 


2 

v 1.?2:;53ri;-02 

3 


. 1370 4E -01 

4 

. 28.138E-01 

o 

1 •'•v:u r ! oo 

1 0 

. 1.5802E+00 

11 


. 1877 IE 100 

12 

.21891E+00 

i. 

oo 

18 

. 79791E+00 

19 


, 91546E+00 

20 

.97959E+00 


'■ Intri-) 

M 

. 31 187E+00 

27 


. 39568E+00 

28 

. 34820E+00 

33 

'll 

:l 

■ooh;-f 

.14 - 

. 11046 E 01 

35 

0 

• 

36 

0 . 

> ■•’..;yEi 02 


.27074E+02 

3 


, 656 7 9 EM 01 

4 

. 46699E +01 


. | 

to 

. 17176E+01 

11 


. 1.5502E + 01 

12 

, 14258E+01 

1 ' 

, ! ’ ■ I- , ! 

1 0 

. 1.0 008 EM 01 . 

1.9 


. 99895E+00 

20 

. 10001E+01 

2 r ‘i 

::i ! oo 

Ms 

. 95585E+00 

27 


.93874E+00 

28 

•92673E+00 

<1 1 

si; to 2 

3 1 

, 9 8 700 El 00 

35 


. 12747E + 01 

36 

. 16415E+01 


MT= 6 

wr= 

• 24435E+01 


1 10B'- . 19572EE008I.)MI\FF = 



2 

, 16180E 02 

3 


. .1 .1 6 4 9 E - 0 1 

4 

. 24729E-01 

9 

, " • ***';'" . * *• > 

1 0 

. 1. 5360E+00 

11 


. .1 8439E + 00 

12 

. 21 A78E+00 

i ■ 


! 

, 79 80.1 EM 00 

.1 9 


, 9 1 5 4 9 E 100 

20 

. 979A7E+00 


■ 1 Mi'. - . . 

..* > 

.3119AEM 00 

2 7 


. i .'MM MOO 

28 

. 34820E+ 00 

r 

■ o !. 

■ 4 - 

, 9741 3E 02 

35 

o 


36 

0 ♦ 


r , ;> 


. 78 34 8E +02 

3 


. 6 “'365EM 0.1 

4 

. 4/312E+ 0 1 


I ' : 4:: v i ! o i 

i 0 

. 1. A343E+01 

1 1 


. 1 16,'-.? EM 01 

12 

. 13439E+01 


■ ■■■:••," lt-1 O'" 1 

1 ■ > 

, 9767 4 E +00 

19 


.98 442E+00 

20 

■ 99453E + 00 


■■I.:-; .1 , 

.V 

9706 4E +00 

2 7 


-'l.; E’iMOO 

28 

. 93985E+00 



• r 

, IvOMEMOl 

35 


1 '.'8 1 >01 

3 6 

, 1. 69 37E' 10 1 


. 13941 1 " 0 1 . 


5 

.59242F-0 1 


. 73569E--0.1 

7 

. 9084 3 E "1 

' * 

1 | i . I f , I . 

1.3 

. 2515.1E EOO 

14 

• 31 1.42E+00 

15 

.3A552E 1-00 

’ ■ 

. 1< ■ :> •! 1 '!• " 

2.1 

. 1.000 OEM 01. 


.98028E+00 

23 

,917013 1 

> 4 

, 800991: I'".) 

29 

• 2861 3E 100 

30 

.24240E+00 

31 

. 18034E1OO 

3 ' ’ 

-V7W6E ■' 1 

37 

0. 

38 

0 . 

39 

0. 

10 ( 

. 

5 

.32042E+01 

6 

. 28239E+01 

7 

.2492 7E. t-01 

8 

.22 1012:.. to 1 

13 

. 13854E+01 

14 

. 12514E+01 

15 

. 117241 i 01 

1 6 

, 1 OH 34 EM 0 1 

21 

. .10 00 OEM- 01. 

. -> , , 

. 99A48E+00 

23 

. 987 1 or:, too 

7 4 

-97843E+00 

29 

. 0825.1 E 100 

30 

. 85 .1 65 E + 00 

31 

. 80873E1 Oo 

33 

.784 HIM 00 

37 

.2 01 5 9 EM 01. 

38 

. 35A39E+01 

39 

• 54983E K) 1 

40 

.25 6 18 EM 02 


.27556E-01 


5 

. 55193E- 01. 

A 

.6927 IE -0.1 

7 

. 8A442EE 0 1 

8 • 1070 OF +00 

13 

. 25049E+00 

14 

. 31 153E+00 

15 

.36 61 OEM 00 

16 . 49R35E 1 00 

21 

. 10000E+01 

'i 'i 

. 98026E+00 

23 

, 9.1.697E+00 

?4 . 80 099 E 100 

29 

. 28558FM00 

30 

. 24.148E+00 

31 

. J.79.13E EOO 

i .. , 96353E-01 

37 

0. 

38 

0. 

39 

0. 

40 0. 

5 

. 31879E-M 0.1 

A 

. 27935E+01 

7 

.2451 OEM- 0.1 

8 .21A06E+01 

13 

. 13299E+01 

14 

. 12016E+01 

15 

.11279EM01 

16 . 10496E 10 1 

21 

. 10000EM 01 

22 

. 1 00 1 2E+01 

23 

. 99547 E EOO 

84 . 9202’tE |O0 

29 

. 89258 E + 00 

30 

. 86015E+00 

3.1 

• 8 157 3 EE 10 0 

3 8 , . '7 1 401: 

37 

.20705E+01 

38 

. 3A755E+01 

39 

.56814E + 0 1 

40 ■ .765 4 1 EM 02 


•40094E- 01 



. 50108E ■ 0 .1 

6 

. 63765E--01 

7 

.0071 4E- 01 8 

. i oi 4 (-00 

13 

.24943E 100 

1. 4 

. 31200E+00 

15 

.36 7. 15 EM 00 .1.6 

, -V-’^O-IE + OO 

21 

.10 000 EM- 01 

':> 9 

. 90032E+00 

23 

. 91700E+00 24 

.80 107 EM 00 

29 

.285 1. OEM 00 

30 

.24065E+00 

31. 

. . 1 . 7809 E 1 00 58 

• , 'o 'f; -o 1 

37 

0 . 

38 

0 . 

37 

0 . ■ 


5 

. 31659E+01 

6 

. 2 750 7E +01. 

7 

.23916EM01 8 

. 20090E-t-0 1 

.1.3 

. 125 1.6 EM 01 

14 

. 1. 1321E+01 

15 

, 1 0659 1:. 1 * 0 .1 

1 002 ir: ioi 

21 

. .1.0000E 10 1. 


. 1. 0072E+01 

23 

. 10056FM 01. .9 4 

1 "045E.+0 l 

‘2Q 

, 9 08 5 9 EM OO 

30 

. 86809E: +00 

3 1 

.021 72 EM 00 v ' 

( 00 

x '? 

. 81 4 441-1 0 1 

b 

.38393EM0I 

39 

.594161 M'M 17 

> J":r I .9-* 


■f 


» 


a > 

M T~ 7 

UT= 

. 26180ET01 

• .1 

MOV . 1846517TOOSUMKFF= 


1894617-01 






f 


9 

. 1. 2661.E- 02 

3 

,~2645E-02 

4 . 2049517-0.1 

5 

.4336917-01 

6 

, 56304E-0.1 

7 

, 7278917 01 

8 

. “3,>9RF -01 

' > > i ■ • < r 

1' 

> i i .v.'f: f oo 

n 

, 17970E100 

12 . 2.1 383E100 

13 

. 24805E100 

.1.4 

. 31269E100 

15 

• 36856FT00 

1 6 

. >r 

i: ■ ■ '’ vb : 

1 <t 

. 798 1 7ET00 

19 

. 1 554E+00 

20 .9797717100 

21 

. 1. 00 OOF. 10 1 

22 

. 9803817100 

23 

,91 701 El 00 

7 1 

.80113F 100 

!./ r'' '-r 


. 51 1. 97ET00 

2 7 

,3959717100 

28 .3479017100 

29 

. 283 99 E 100 

30 

.2389017100 

31 

. 175 91 El 00 

37 

,94 20 IE -01. 

, : ; M.3nc ■ -oi 

•' i 


■ .6766317-02 

35 

0, 

36 0. 

37 

0. 

38 

0. 

39 

0. 

40 

:) . 

: . ; n'oi:-; -.>■> 

9 

. 299.1. 3ET02 

3 

,6927317101 

4 .4785817101 

5 

. 3 1 1.3317 1 0.1. 

6 

.2671.31710:1. 

7 

.2291017101 

8 

. 19732ET01 

!. :I. E ’■ l 

.10 

• 1506.217101 

1 .1. 

. 13390E101 

12 . 122.1.7ET01 

13 

..1.136117101 

.1.4 

. 10309ET01 

15 

. 97637ET00 

16 

.9328717100 

:! ”■ ■ ;t;-? oo 

111 

• ? 1 < 1 217100 

19 

. 96 478 ET 00 

20 . 98713ET00 

21 

.1000017101 

n 

. 1014517101 

23 

. 1. 0 .1 7 7 E 1 0 1 

24 

. 1.020? 1710! 

2‘ - : I'VH '/l-'-fOt 

.V, 

.98685ET00 

27 

, 96 77 4 E 1 00 

28' .95176ET00 

29 

. 90966E 1 00 

30 

. 87212E100 

31 

, 8 232 4 FI 00 

r n 

,8008117100 

3 ■ 1 1 : - 

4- ■ /■ ri •: 

34 

, i0324i7ioi 

35 

. 1 35 06 E 101 

36 ♦ 17582E10.1. 

37 

.2236317101 

38 

. 40265E101. 

39 

. 626)421’ 1 0 1 

10 

, 29532ET02 

nr- 8 

WT 

= . 27925ET01 

HDD- . 17436I7T00SIJMKFF= 

.51 364E-01 







:i. ■ f 

o 

. 8768417-03 

3 

.6564117-02 

4 . 15280E-01 

5 

.34306E-01 

6 

. 45920E-01 

7 

.6138317-0:1 

8 

.8160:.H7 -01 

o } r% "•.)/, |;r j .-,r, 

1.0 

. .1.3 771 17100 

.1.1 

. 1 7261ET00 

12 . 20949E100 

13 

. 24607E100 

14 

. 31358ET00 

15 

. 37042ET00 

1 6 

, 5003 1171 00 

1 ••• :v>m ! .U U'O 

18 

. 79834ET00 

.1.9 

,91.56517100 

20 . 97983E1 00 

21 

. .1000017101 

22 

. 98037E100 

23 

.9169117100 

24 

.not. one too 

; v : ■ :>o\ r:'h roo 

26 

. 51 .1. 81 17100 

27 

. 39589E100 

28 .3471SET00 

29 

.28135E100 

30 

. 23495E100 

31 

. 171.1217100 

32 

.9093717-01 

T3 ''372! II. 01 

S 1 r 

:m ■ 

. 3097317-02 

35 

0. 

36 0. 

37 

0. 

38 

0. 

39 

0 . 

40 

0 . 

1 . A7O74! :; i07 

9 

. 31508ET02 

3 

. 70720E101 

4 . 47803E101 

5 

» 29861ET01 

6 

.25132ET01 

7 

• 21089FT0 l 

8 

. .1 7750E101 

? 5 7! 2 8 13! 01 

1.0 

. 13006ET01 

11 

. 1 1.400E101 

12 . 10342E101 

13 

• 96184E100 

14 

. 881 76E100 

15 

.84607E100 

16 

.8329917100 

1 7 V'.33l- f 00 

1.8 

. 89998ET00 

19 

. 94.L50ET00 

20 • 97791E100 

21 

. 10000E101 


. 10217ET0.1 

23 

. 10 288 El 01 

24 

. 10355E101 

3'* ; !••!. 3/>r-J01. 

26 

. 9 97 52 E ( 00 

27 

, 97403E100 

28 . 95281E100 

29 

. 90233E100 

30 

. 86060E100 

31 

.81002E100 

32 

.7926917 TOO 

->: \ • ; | :; { > 
.1 | VF .5 f 

31 

. 1055117101 

35 

, 139 13 El 01 

36 . 1823.1.ET0.1. 

37 

.2330217101 

38 

. 42288E101 

39 

.660101710,1 

10 

.3125317102 

MT- 9 

WT 

= . 29671ET01 

HOD = . 1.6702E100SUMKFF- 

. 48080E -0 .1 







i 

9 

. 4554 1. E— 03 

3 

.35265E-02 

4 . 88070E-02 

5 

.2164917 0 1 

6 

.3059817-01 

7 

. 4355817 -0 1. 

8 

,6216917 -0,1 

0! 

I 0 

. 12107E100 

11 

. 1 6077ET00 

.1.2 . 20291 E100 

13 

.243691 100 

14 

. 31540E100 

15 

. 37382F 1 00 

1 . 6 

.5017517100 

i. • yy ko 

.1.8 

.7987917100 

1.9 

.9159217100 

20 , 97983E100 

21 

. .1.000017101 

,.j 

. 98026E100 

23 

. V.1.672ET00 

21 

. 8009017 100 

• f 00 

.'■■A 

.5U63ET00 

27 

, 3955217100 

28 .34521E100 

29 

. 27656 E 100 

30 

. 22 84 2E TOO 

31 

. .1.6 4 08 El 00 

32 

.89 0 4 317 -01 

, • 7 -7 ( '. '■ ,! ■: -Oj 

4 1 

; 1 1 • 

- . 1 3270E-02 

35 

- , 44577E-02 

36 0. 

37 

0 . 

38 

0. 

39 

0 . 

40 

0 , 

i . ' 'r'." i :i :• 

V> 

.3269SE+02 

3 

. 70 6 ft 4 El 01 

4 . 46389ET01 

5 

. 27245E10.1 

6 

.22199ET01 

7 

. 17912E101 

8 

. 14430F10I 

,■ - | MU 

. 1 . o 

. 97742ET00 

11. 

, 83833 El 00 

.1.2 . 76028ET00 

1.3 

. 7161 9ET00 

1.4 

. 68451E100 

15 

.6 8363 FI 00 

1 6 

.7 239017100 

1 ' " ! >or ; Or' 

111 

, 86 76 9E TOO 

19 

. 9260617100 

20 .97127ET00 

21 

. 10000E101 

22 

. 10246E101 

23 

. 10323 El 01 

74 

, 1038517101 

3'. ■ -■■■■■-:■: Of- ! •• ; 


: 9 8 <? 5 2 E l 00 

27 

, 95475E100 

28 . 92350E100 

29 

, 8 A 05 OE 100 

30 

. 81 434E100 

31 

. 76 389F100 

.17 

,■75527171 00 

'.13 •••; • » ? :u }<•«> 

VI 

. 104 77E10.1. 

35 

, 1 398017101 

36 . 18546E101 

37 

. 23S80ET01 

38 

, 43771ET01 

39 

, 68585E1 01 

40 

. 326iO IE 102 




WT= • 31 416E+01 HOB- . 16329E+00SUMKFF= . 4256 It: -01. 




.584 "OF! -05 

3 

. 1 0905E 

•1 

6 .1 054 E~ 03 


7 ('). 

1 0 

. 8S1.58E -01 

1 1 

. 14990E+00 

12 

.21019F1 00 

13 


1 0 

.?9?:v?Ei00 

1° 

. 916.1.4F+00 

20 

. 9 7900 E +00 

21 


6 

. 0 1 1 .",•>+ ! 00 

27 

• 39255E I 00 

28 

. 33693E+00 

29 

• •'» 1 

34 

.17989E-01 

35 

. S0616E-02 

36 

.37814E-02 

37 

o\: ny? 

2 

• . 3327 2 E +02 

3 

.69053E+01 

4 

• 43652E+01 

5 

?p '■ () n 

1 0 

.59241E+00 

1.1 

• 54258E+00 

12 

. 55458E+00 

13 


If? 

.871 55 E+ 00 

19 

> 929.1. 3 Ei 00 

20 

.97231E+00 

21 

oo 

0 

< '’4459E+00 

27 

. 86909EI 00 

28 

.81 336E+00 

29 

f>l\ * Oo 
1 r • : •. •"> 

; a 

.10P39E+01 

35 

. 1 1426FF01 

36 

. 1S992E+01 

37 


26038F -0? 

6 

. 5 0 1 2 3 E - 0 2 

7 

. 999541 1! 

. ?,>74.'.|;: 0 1 

259361. 1 00 

14 

. 33192E+00 

1 5 

. 38624E 100 16 

• 3071 3E + 0O 

1000 OF! + 01 


. 98007E+00 

23 

.91 6601+00 74 

800:,.'+; * oo 

26015E+00 

30 

. 20804E+00 

31 

.14546+100 37 

. 8324.1 1 o 1 

18259E -02 

38 

. 37005E-03 

39 

. 5121 3E -04 40 

. 2 4 407F-03 

23353E+01 

6 

. 17957E+01 

7 

.1339 31)0! 

. 9,204 1 + 100 

58024E+00 

14 

. 63027E+00 

15 

» 66558E t 00 16 

. 73000E+OO 

lOOOOE+Ol 

22 

• 10199E+01 

23 

. 10236E + 0 1 24 

. 1019(i|; in | 

731.1 6F+00 

30 

. 68804E+00 

31 

. 661 39F 4 oo 32 

. 7 01 8 2E +00 

24387E+01 

38 

. 44673E+01 

39 

.700501 I'M 40 

.33353E102 


WT- . 33161E+01 HOB- . 161 1 6E+00SUMKF F~ 




. 45178E-03 

3 

. 34576 E- 02 

4 

.83580E-02 

■ 

•12o;v5E--01 10 

. 1 1 439E+00 

11 

. 153 78E+00 

12 

. 19609E+00 

1 

' "-'F F '">0 12 

. 7 9999E+00 

19 

. 9 1. 657E + 00 

20 

. 97997E+00 


■ ' v+oor.. : oo 26 

. 30S38E+00 

27 

.3891 7E+00 

28 

.33963E+00 

,«v 

551 7 4F-07 24 

- . 10609E01 

35 

0 ■ 

36 

0 . 

1 

' ’51 4E F02 2 

. 33891E+02 

3 

. 72946F+01 

4 

.47637E+01 


' : 5 12 +’ (01 I 0 

. 9 1 3 6 4 E 4 0 0 

1.1 

. 76730E+00 

. 1.2 

.68709E+00 


, 23383F ! 'H> .18 

. 95557E+00 

19 

. 994 1.6E + 00 

20 

.10065E+01 

*■ ' ■’ 

87 84 9+ +00 26 

. 79506E+00 

27 

. 75049E+00 

28 

.73014E+00 

7 o 

* i 

> A [ 1..' 1 1 > 

V--OFIO? 

. 87478E+00 

35 

. 12637E+01. 

36 

. 17471E+0.1. 


- » 19066E-01 


5 

. 201 75E-0.1 

6 

. 28388E-01. 

7 

.4034 OF O' 

(!{ 

» 5 /6m'E -0 l 

13 

. 23729E+00 

14 

. 31001E + 00 

15 

. 37059F 100 

16 

♦ oo i '£■} oo 

21 

. 1. OOOOE + O.l 

'■> 

. 97982E+00 

23 

.916101! 100 

2 1 

. 7°? Vf F: f oo 

29 

. 27716E+00 

30 

. 23445E.+00 

31 

. 17122EF00 

3 : 

* -701 'V.L - 0.1 

37 

0 . 

38 

0 . 

39 

0 . 

40 < 

) 

5 

• 2761 1E+01 

6 

. 22305E+01 

7 

. 17778FM01. 

8 

• L - 10 B/F: f 01 

13 

. 644 73E+00 

14 

. 62455E+00 

.15 

. 63928F+00 

16 

. 1 or* 

21 

. 10000E+0 1 

2 2 

. 98282E+00 

23 

. 94915E+00 

24 

. ^0A77F. f 00 

29 

. 691 48E+00 

30 

. 65550E+00 

31 

.601 15F ! 00 

' I , , 

+ 7 'o 7 *ir- j-oo 

37 

. 23079E+01 

30 

. 43884E+01 

39 

.69740+ O' | 

40 

. 7771 0 1 ;: Hv;> 


MT-- 12 


WT= > 34907E+01 


HOI'I- . 1S44OE+0OSUMKFF-- . .1 8799E- 0 I 


. 8 7 5 6 3 F! ■ 0 3 



; O.'o;.! 


10 

. 1 381 RE + 00 



: f 00 

1 8 

. '>97 4AE+00 


■■'481 

r '■ > ,y 'i 

0 A 

. 509 42F+00 

77 
1 ! • 

■ . ■:i| 

. .»> 


) . 

l 

| , 



- 15494E+02 

Q 

1 59 171 

f i 

i n 

. 1 3 -15 6 El- 01 

1 '• 

' ’ 5 1! 

: o 

X 8 

. 8 1365E f 00 

i: 


7 00 


. 8 ? 884F: + 00 


' ■.. . • .;■> 

.4 A 

. 10106 E+ 01 


•"> 

- 63575E- 02 

4 

. 1 51 26E -01 

11 

. 1 7350E+00 

12 

.21056E+00 

19 

. 91550E+00 

20 

.97977E+00 

27 

. 39048E+00 

28 

.341 76E+00 

35 

0 . 

36 

0 . 

3 

. 78826E+01 

4 

•52891E+01 

11 

. 1 1627E+01 

12 

. 10404E+0.1. 

19 

* 924,691!;+ 00 

20 

. 97657E + 00 

27 

. 86039E+00 

28 

.84739E+00 

35 

. i 390 !f FO 1 

36 

. 18902E+01 


5 

. 339 4 71; -ot 

6 

.45508F -01 

1 3 

. 24 70 OF 100 

14 

. 3.1 377E+00 

21 

. 10000 E +01 

9 ■:> 

. 98011E+00 

2 9 

.2814 51!! 100 

30 

. 24003E + 00 

37 

0. 

38 

0. 

5 

. 32567E+01 

6 

. 2720 8:. 101 

13 

. 95 4 5 OF! + 00 

1 1 

. 851 6 OFF 00 

21 

. toooomoi. 

7 2 

. 10085E + 0.1 

29 

. 81767E-! 00 

30 

. 78406E+00 

37 

.24673+ 101 

38 

. 16.22 41!. + 01 


7 

.60982+ v; 1 

" 

, 8 1 .5 ■ i ) 1 

15 

. 36773+ t 00 

1 .', 

,4984.5+ 100 

23 

. 9 1 6,57 E *80 

' ' 1 

■ 80 00°t! 100 

31 

. 1 769 4E 1 0" 

T ’ 

.8404 IF -01 

39 

0, 

4 0 


7 

. 22622+ ! 0.1. 

7 

, .1.88,58+ 1 0 1 

15 

.79437+1 no 

1 

, -" 1 4 57 + i O ' 

23 

.993361 5,. . 

I 

9 " i 7 5! 1 

31 

. 727621 * 00 

7* 

.69 295+1 00 

39 

. 730981- 4 m | 

V. 

'■ , ! r" ( O'' 


I 


0 


MI 1.3 WT= . 36652E+01 HOP- . 1. 4171E +OOSUMKEF- 




. .1 3072E-02 

3 : 

97123E-02 

4 

.211 7 4 E - 0 .1. 

riO'-: ! ■■ 

) 1 n 

. 152 19 El- 00 

.1 1 

1 85201:: loo 

12 

. 21977E+00 

■ - hi 

■i ! H 

r i o o 

1 9 

9 1.473!;: (•()() 

20 

. 97946E+00 


•' A 

.3 1 200 E-10 0 

27 

393 70 Ei 00 

28 

. 34525E4 00 

'•"■I- '■ 

| 'A , 

. 

35 0. 


36 

0. 

•; (.:• r<- 


< 38 7 1 OE 1 02 

3 , 

99301 E l 01 

4 

. 60552E+01 

: 'l ;; ; .'V 

t 1. o 

■ .1 0 >''31.101 

11 

1 A 283 E + 01 

1.2 

. 14787E+01 

1 >\i\' I 

! r; 

.9.1 983E + 00 

19 . 

930.521::+ 00 

20 

. 9681 1E+00 

i : : o 

1 :6 

. 1 00071:: t o i 

27 . 

9762 IE +00 

28 

. 9674RE + 00 

/. 

•i or*. 4 


. 1 1 067E + 0 .1. 

35 . 

1511 0E+01 

36 

, 20336 E+ 01 


M1 14 WT= . 38397E+01 HOP- . 13966E+00SUMKFF- 



. 16643E-02 

3 

12167E-01 

4 

•25783E-01 

10 

. 1. 6052E+00 

11 

1. 9207E+00 

12 

. 22489E+00 

• ■ * ; ■;> !. n 

. 79906E+00 

19 

9 156> IE +00 

20 

. 97972E+00 

26 

. 5.1. 3 47 E fOO 

27 

39560E+00 

28 

. 34669E+00 

, >] -■■)[ > ■) 

:>. 

35 0 


36 

0. 

■■A ) ! . ; « ;• '> '» 

, 39331E+02 

3 

91756E+01 

4 

. 63958E+01 

io 

.2:16211- +01 

1.1 

19373E+01 

12 

. 1 7744E+01 

> ' •’ 1 . ! ■ ' 1 is 

. 1 1008E+01 

19 

10369E+01 

20 

. 101 13E+01 

i i.o;:: ! - > i : v> 

. 10049E + 0.1 

27 

98770E+00 

28 

.98187E+00 

. • i ' ;m 

. U046E + 01 

35 

15025E+01 

36 

. 20236E+01 


HT' = 15 WT = .40143E+01 HOP- . 15446E+00SUMKFF~ 


2 

, 19 7 

59E 

- 02 

3 

1. 1278 E 

-01 

4 

. 29205E-01 

1.0 

. 163 

7 E 

+ 00 

1 1 

1.91021: 

100 

12 

. 22558E 1 00 

1 9 

'7 Q Q P: J jo p 

+00 

1 9 ■ 

. 9 .1 6 1 5E 

100 

20 

. 98002E+00 

’•"> 

. • ; i :> 

70 E 

+ 00 

27 • 

. 39476E 

+ 00 

28 

.34522E+00 

‘i 1 < 

> . 



35 0, 



36 

0. 


.357 

82E 

+ 02 

3 « 

•85545E 

+ 01 

4 

.60491E+01 

! « •' 

.219 

7 0 E 

+ 01 

1 1 

, 19825E 

+ 01 

12 

. 18236E+01 

1 S 

. 1 1 6 

49E 

+ 01 

1 9 

, 1091 7E 

4 01 

20 

, 10441E + 0.1. 


, 9 1. 5 

3 6 i: 

+ 00 

27 

■ 90006E 

+ 00 

28 

, 89469E+00 

1 1 

, 101 

69E 

+01 

35 

. 13 80 2E 

10.1 

36 

. 18537E+0 1 


.24493E--01 


5 

. 4 4 705 E - 01 

6 

.58088E -01 

•7 

. 7 5 .1. 7 2 E: -01 

8 

'-V.532E- 0.1 

13 

. 25398E 100 

1.4 

. 31780E4 00 

15 

. 371 59+ 100 

I 4 

707 1. 7i: : ' •! 00 

21 

.100001: 101 

') 

.98060E 100 

23 

.917181:+ 00 

.7 1 

801 111:. +00 

29 

. 2 872 3 E 4- 00 

70 

• 24862E +00 

31 

.18967F 400 

<7 

99,5 19+ -3) 4. 

37 

0. 

38 

0. 

39 

0. 

10 0 


5 

. 309 5 4 +.'401 

6 

. 33270E+01 

7 

.203951: 40 1 

8 

2 13381: 401 

13 

. 13679E + 01 

14 

. 12253E+01 

.15 

. 1.1428+ ‘ 0 1 

16 

I0516E4 04 

21 

. 100001: 101 


. 103 HE +01 

23 

. 103961: 4 0 I 

24 

1013 vi: 4 0 4 

29 

. 9461 OE +00 

30 

.91670E+00 

31 

.86064+ 400 

7:7 

81 683 E 4 00 

37 

.2651 IE 101. 

38 

. 49722E+01 

39 

.78795+401 

40 

38157E+02 


.26756E--01 


5 

. 5240.1E-01 

6 

. 66772E-01 

7 

. 84592!:. 0 l 

<:) 

. .1 062.5+: ! 00 

13 

. 25763E4 00 

14 

.31965E+00 

15 

•37352EI00 

1 6 

.50115+ 100 

21 

. 1 OOOOE+O 1 

22 

. 98060E+00 

23 

.91736E+00 

2 4 

.80176E 400 

29 

. 29028E+00 

30 

. 25447E+00 

31 

. 19949E+00 

5 2 

, ! 1072+100 

37 

0. 

38 

0. 

39 

0. 

40 ( 


>.) 

. 42391E4 01 

6 

. 36710E+01 

7 

.31821E+ 01 

8 

,2/72’AE i 0 1. 

13 

. 1650BE+01 

14 

. 14076E+01 

15 

. 13922E+01 

16 

. 12808E4 01 

21 

. 1 OOOOE+O 1 

'>':) 

. 10091E+01 

23 

.10118+101 

24 

. 1 02 23 E 4-0 1 

29 

.97028E+00 

30 

. 94867E+00 

31 

.89882+ +OO 

32 

» 8307'.- +: 400 

37 

. 26)422+".+ 01 

38 

. 49785E + 01 

39 

.79132E4 01 

40 

.3 856 IE +02 


- . 3471BE-0 1 

5 .5733.1+ 01 6 .71991E-01 7 .B9053E--01 0 < I 1 1 1 9 1" 4 0 0 

13 . 25 7 2 7E -100 11 .31B02E+00 15 .371771 100 l A 50271E+00 

21 . 1. 0000 E 101 22 . 98027E+00 23 .91692+400 24 M40 .I.08E 4 00 

29 . 28905EI 00 30 .25430E+00 31 .20023E+00 33 , 1 0982+. 1 «0 

37 0. 39 0. 39 0. 10 0. 

5 .4103 IE 10 l 6 . 3588 IE 10 1 7 .31427E+01 0 , 27 663+: I 0 I 

13 . 1700 VE I'M. 14 , 15353E10.I. 15 . 14372E4 01 1.6 . 1 52 I OE 40 I. 

21 , .1000 OE 101 22 , 9744 IE 100 23 .95086+: I 00 31 . T I 

39 , 0974 IE ! 00 30 .87094E+00 31 .82784+400 52 . :, 82i:l7E f-o 

37 .2415 1 E I 0 1 38 .45330E+01 39 .71 9131. 1 01. 10 . ’ +• M' '07 


mt= 16 

WT= 

.41888E+01 

HOB . 17601 E+OOSUMKFF= 


96244E-01 







;> 

. 72.1. -12E -02 

:j 

. i::;799E-oi 

4 

. 3 1.5 1.9 E" 01. 

5 

. 6036417-0 1. 

/, 

. 750 4917-01 

7 

. 92745E "01 

8 .113671:100 

i 3775 e : oo 

10 

. 16472E100 

1 1. 

. 1940317100 

12 

. 22466E100 

1.3 

. 25557E+00 

14 

. 31539E100 

15 

. 369051: 100 

1.6 . 5006 7171 00 

i ■ 5^75'ir )~oo 

18 

. 79906E+ 0(> 

1.9 

• 9.1602F 100 

20 

. 97996E 100 

21 

. lOOOOEtOl 

22 

. 9B005E+00 

23 

. 9166 1 E! 00 

• 1 .800441 100 

, 60O56E K ' 0 

26 

.511 50E100 

27 

. 3931517100 

28 

. 34380E 1 00 

29 

. 28720E f 00 

30 

. 251 19E100 

31 

. 1 9470F 1 00 

57 .102431:100 

. 1 0822E-01 

A 1 A , 

3-1 0 


35 

0. 

36 

0. 

37 

0. 

38 

0. 

39 

0. 

40 0. 

; 54 417102 

-;> 

.31 602E+02 

3 

. 7 71.56E-101 

4 

.55162E+01 

5 

. 30033E 101 

6 

. 33478E101 

7 

. 295 2 OE 10 1. 

8 .26.1.571-10:1 

: 33 1 ?'/f : 01 

.1.0* 

. 20994E+01 

i. i 

. 19021E+01 

12 

. 17536E+01 

13 

. .16375E+01 

14 

. 14782E+01 

15 

. .1 38281710 1 

1.6 .176901101. 

1.7 , J.2099E101. 

18 

. 11393E101 

19 

.10866EK11 

20 

. 10434E+01 

21 

.lOOOOEtOl 

22 

• 96496E+00 

23 

. 9298.1. E 100 

.74 . 90242E 100 

.871 06 FI 00 

26 

. 861 1 7Et00 

27 

. 84363E+ 00 

28 

. 83799E100 

29 

. 82870E 100 

30 

.81121E+00 

31 

.7699317100 

32 .730971100 

. 7 70 4 OF 100 

34 

• 95651.E+00 

35 

. 12B20E-101 

36 

.17022E+01 

37 

. 21990E101 

38 

• 40669E 101 

39 

. 6407 OE 101 

40 .307821102 


MT= 17 UT= . 43A33E+01 MOB= . 194B5EtOOSUMKFF== -.92078E-01 



2 

. 23866E-02 

3 

. 16874E--01 

4 

. 331 15E-01 

5 

. 62396E -01 

6 

.77 077E-01 

7 

.9464517 01 

a 

.115291" 100 

1 3893E+00 

. 1.0 

. 16533E+00 

11 

. 19400E400 

12 

. 22401E100 

13 

. 25440E100 

14 

. 31354E100 

15 

. 367081 K'O 

16 

. 4991.1.1100 

77-:f 151)00 

18 

• 79852E+00 

19 

• 91576E400 

20 

• 97982E+00 

21 

. lOOOOEtOl 

22 

. 98006£t00 

23 

.91.6631 10O 

>4 

. 800381 1 00 

f 00 

26 

.51132E+00 

27 

. 39317E100 

28 

. 34455E 100 

29 

.287891 100 

30 

. 25.1 53 El 00 

31 

. 1971.91 'OO 


, 1. 1. 1 1 217100 

a o580E. -02 

34 

~ . 3.1. 938E-01 

35 

0. 

36 

0. 

37 

0. 

38 

0. 

39 

0. 

40 

0. 

6 1 232E102 

2 

. 28678E+02 

3 

• 71061E+01 

4 

. 51 185E101 

5 

.3567717101 

6 

. 31538E101 

7. 

.279311710 l 

8 

. 2 4 85 3E 10 1. 

2 93551710 1. 

10 

. 200B7E+01 

11 

. 1 8247E 101 

12 

. 16849Et01 

13 

. 15746E10.1 

14 

. 1. 421 7E101 

15 

. 13294F401. 

16 

. .12181 El 0 l 

1 1. 6912171 01 

1.8 

. 11046E101. 

19 

. 10654E+01 

20 

. 10333E101 

21 

.lOOOOEtOl 

22 

, 97034E100 

23 

.937707100 

2 4 

.9084017100 

8 725 8 E 100 

26 

. 85998E100 

27 

. 34087E+00 

28 

. 83467E 100 

29 

.82246Et00 

30 

. 80420E100 

31 

.7661817100 

32 

. 72626 E 1 00 

73436EIO0 

34 

.8707217100 

35 

. 11893E+01 

36 

. 15853Et01 

37 

.2042017101 

38 

. 37422E 101 

39 

. 586321710 1. 

40 

. 27902E 102 


'.•'■6 K'2 

MT ~ 18 WT= . 45379E+01 HOB* . 2081 7E100SUMKFF* -.740B8E-01 



2 

♦ 1M982E-02 

3 

, 1757417-01, 

4 

.34 118E -01 

5 

. 6364 c>F -01. 

6 

. 78324E-01 

7 

.9581 T 

• i\ I 

<! . 1 .1 6 30F 100 

'0! iO() 

10 

. 16578E100 

11 

.1940817100 

.1.2 

.22372E100 

13 

.253 00 E 100 

14 

. 31254E100 

1.5 

.365991 

■HH) 

16 .49821 1-100 

1 T! 00 

to 

. 79S23E100 

1.9 

. 91 559E100 

20 

. 97975E 100 

21 

.lOOOOEtOl 

92 

. 980.1 2E100 

23 

.916 7/1. 

1 Ot< 

' 1 .80053F-100 

■ ,v: T 1 00 

2 6 

, 51 1 35E100 

27 

. 3935317100 

28 

. 34570F 100 

29 

.2878317100 

30 

. 24887E100 

31 

.1910 11- 


'. , . 1 0459i:: 100 

1 '.I:. -01 

34 

• . 25156E-01 

35 

0, 

36 

0. 

37 

0. 

38 

0. 

39 

0. 


40 0. 

/• } } n .* 

'■> 

. 2691 9 E l 02 

3 

.6726117101 

4 

. 48646E-101 

rr 

.341.0217101 

6 

. 302.1 2 Et 01 

/ 

• 248 1 61: 

1 V 1 

■ ; .33:9-1 11 -101 

4 1 ! 0 1 

1. 0 

. 19 3 93 El 01 

1 1 . 

. 17639E101 

12 

.16299E101 

13 

. .1. 5237E 101 

14 

. 13757E101 

1 5 

. 1 286 1 E 

10.1. 

1.6 , 1 1 .’00 El 01 

75 1 1 0 1 

1 8 

. .1.0 7 74 EM 01. 

19 

. 1 0 475E101 

20 

. 10246E101 

21 

. 10000E 101. 

9 9 

. 9781. 7E tOO 

23 

.952461: 

1 00 

, ,, 99 * /nr f 00 

\ jj- -j i’i 

26 

. 38050E100 

27 

. 8620 4171 00 

28 

. 85534E100 

29 

. 83 74 4E 100 

30 

. 81 478E 100 

31 

. 7716'M 

l 

. “,19651 too 

•4 ” i: : ; •* •><> 

34 

. S9556EtOO 

35 

. 1 1 90 T 101 

36 

. .1.5588E 10.1 

37 

.198571" 101. 

38 

.3576617101 

39 

.5 56 2 IF 

f 0 1 

10 • 96 1 102 




r 


H I' :: 19 WT= . 47124E+01 IIOFi . 21591E+00SUMKFF“ 


1 • 

> ( 

o 

.25296E-02 

3 

. 17753E -01 

4 

• 34376E-01 


l 7 9 ~’2Z (-00 

10 

. 1 6567E+00 

11 

. 19383E+00 

12 

♦ 22335E+00 

:l ■' 

' : 1 46E+00 

18 

. 79802E+00 

19 

. 91548E+00 

20 

. 97967E+00 

') 

4011 PE (00 

26 

.511 47E+00 

27 

. 39420E + 00 

28 

. 34688E+00 

41 ■ 
1 

. 1 'M 

34 

- . 1 9919E-0-1 

35 

0. 

36 

0, 

' ■ OO'-F i -■ 

2 

.25982E+02 

3 

. 65160E+01 

4 

• 4721 1E+01 

■ :> 

‘1 043F 101 

10 

. 18957E+01 

11 

. 17254E+01 

12 

. 15951E +01 

I 7 

: ‘ 0 17E ! 0 1 

1 0 

. 1 0603E+01 

19 

. 1 076 OF + 01 

20 

.. 101B3E + 01 

25 

. • 1 450F +00 

26 

•90017E+00 

27 

. 88287E+00 

28 

.87507E+00 

33 
4 1. 

■Q n,\F)0-> 
•; " .r l oo 

1 

. 9 1882E+00 

35 

. U995E+01 

36 

. 15528E+01 


-.46652E-01 


5 

• 63940E-01 

6 

. 78593E-01 

7 

. 96036E- 01 

8 

. 11644E+00 

13 

. 25331E+00 

14 

. 3.1 192E+00 

15 

.36536E+00 

1 A 

. 4 9 76 9E +00 

21 

< 10000E+01 

':> 

. 98015E+00 

23 

.91686E+00 

24 

. 80067E+00 

29 

.2B739E+00 

30 

. 24632E+00 

31 

. 1863.1 E + 00 

32 

- lOO/VE +00 

37 

0. 

38 

0. 

39 

0. 

40 

0 . 

5 

•33180E+01 

6 

• 29425E+01 

7 

.26143E+01. 

a 

. 23334E+01 

13 

. 14917E+01 

14 

. 13472E +01 

15 

. 12596E+01 

16 

. .1 1539+: I 01 

21 

. 10000E+01 

o r> 

. 98393E+00 

23 

. 96389F (00 

> i 

. 1 4 I'.’! (00 

29 

♦85205E+00 

30 

• 82621E +00 

31 

• 78512E+00 

32 

, 7546 3E +00 

37 

. 19634E+01 

38 

. 34957E +01 

39 

.54090EF0 1 

40 

. ( or* 



APPENDIX B 

DERIVATION OF Q, 
k 





APPENDIX B 


Derivation of Q, . 

k*A 

The kernel functions Q . were developed originally in Reference [ 14] 

k> J 

and are reproduced as follows: 


X k,j 

Vj 






K, (k j)-K (k j) for j = 1 

i o, i o, 

K ? (k,j-1) - K (k j-1) for j = 2, 4, 6,. <, ,,k -1 
[^(k.j) + K 3 (k,j-2)] - [^(k j) + K^k j-2)] 


for j = 3, 5, 7 k f -2 


= K 3 (k,j-2) - K 3 (k Q j-2) for j = k f 


where 


K 1 (k» j ) 26 ^ ("3Vj v j+2 ) " — ^2 ~ u j( £n |u^ | - 1) 
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2 36 j 


K ->(k» j) (-v. - 3v. +9 ) 1 -+ u.^ (to|u.^ 9 | - 1) 
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3 V *' ,J / 26, v v j j+2' ‘ “j+2 v ~“> “j+2 1 


36 . 
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and 


5 j ‘ - h 






APPENDIX C 


r ^k+1/2 . _ 

Derivation of ^ for Equations (37) 

J 

Equations (37) are linearized to form the following system of equations 
which are solved for P^ to P^ _j, P^ to P^ _ 1 » c |, and C q . 

K A o o f f 


^k+1/2 
. V -— Ac! + 
3cj. I 


k -1 
o 

I 

j =k A 


n 


Bill. 


k+1/2 


3P. 

3 


k ff 1 


n 


AP n+1 + 

3 


Z 


r Dtp 


k+1/2 


j=k +1' 
J o 


3P. 

3 


4P „« + ic . . 

J o 


’•'k+1/2 


k = k. ... k ~1 
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where 
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= 0 


for k = k. to k = k -1 
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